{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3d0ef82b",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Replication file for 'The Intertemporal Keynesian Cross'\n",
    "### SHIW appendix\n",
    "\n",
    "Adrien Auclert, Matt Rognlie, Ludwig Straub\n",
    "\n",
    "April 2024\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd671e10-8582-4196-868c-5f02f00c2932",
   "metadata": {},
   "source": [
    "This notebook contains:\n",
    "- Analysis of SHIW data - data for Figure 1 and appendix Figure C.1(a)\n",
    "- Verification of FOSD property for distribution of MPCs weighted by date-0 income for Figure C.1(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "461903cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "import calibration\n",
    "import models_heterogeneous\n",
    "from sec34_plots import set_texfig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc9625ed",
   "metadata": {},
   "source": [
    "Import 2010 and 2016 SHIW data, and keep only incomes above 0:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cfc25513",
   "metadata": {},
   "outputs": [],
   "source": [
    "s2010= pd.read_csv('data/shiw2010.csv')\n",
    "s2016= pd.read_csv('data/shiw2016.csv')\n",
    "s2010 = s2010[s2010['y']>0] \n",
    "s2016 = s2016[s2016['y']>0] "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "037807b8",
   "metadata": {},
   "source": [
    "We have identifiers, sampling weights, incomes, and survey-based MPCs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "972ea733",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nquest</th>\n",
       "      <th>pesofit</th>\n",
       "      <th>y</th>\n",
       "      <th>riscons2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>173</td>\n",
       "      <td>0.259294</td>\n",
       "      <td>33267.445</td>\n",
       "      <td>60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>375</td>\n",
       "      <td>0.771570</td>\n",
       "      <td>14580.000</td>\n",
       "      <td>60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>465</td>\n",
       "      <td>0.308497</td>\n",
       "      <td>57819.023</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>629</td>\n",
       "      <td>0.202026</td>\n",
       "      <td>37564.723</td>\n",
       "      <td>40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>632</td>\n",
       "      <td>0.164572</td>\n",
       "      <td>47306.563</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   nquest   pesofit          y  riscons2\n",
       "0     173  0.259294  33267.445        60\n",
       "1     375  0.771570  14580.000        60\n",
       "2     465  0.308497  57819.023         0\n",
       "3     629  0.202026  37564.723        40\n",
       "4     632  0.164572  47306.563         0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s2010.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dfd5d03b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Number of rows, SHIW 2020 =  7936  SHIW 2016 =  7367\n"
     ]
    }
   ],
   "source": [
    "print(' Number of rows, SHIW 2020 = ', len(s2010), ' SHIW 2016 = ', len(s2016))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0fe5456",
   "metadata": {},
   "source": [
    "Generate an income-weighted MPC CDF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "276d5b95",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_cdf(data):\n",
    "    stats = data.copy()\n",
    "    stats.rename(columns = {'riscons2': 'MPC'}, inplace=True)\n",
    "    stats['yweight']=stats['pesofit']*stats['y'] # Income weight\n",
    "    stats = stats.groupby('MPC')['yweight'].agg(['sum']).rename(columns = {'sum': 'weight'})\n",
    "    stats['freq_weight'] = stats['weight'] / sum(stats['weight'])\n",
    "    stats['cdf_weight'] = stats['freq_weight'].cumsum()\n",
    "    stats.reset_index(inplace=True)\n",
    "    stats.loc[-1] = [-1, 0, 0, 0] # Add first data point for plotting\n",
    "    return stats.sort_values(by='MPC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "35f44f23",
   "metadata": {},
   "outputs": [],
   "source": [
    "stats2010 = gen_cdf(s2010)\n",
    "stats2016 = gen_cdf(s2016)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0f7179d",
   "metadata": {},
   "source": [
    "# Figure C1(a): income-weighted MPC distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "922b9d2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEaCAYAAADzDTuZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAkJElEQVR4nO3dzW8j550n8O/PMRoYt1+qJXfQPmxbTQ2S8SGAhlKM2VMDDmXMIQMEG6ob6JOxmKZmb3PIku6dP6BDzRwGyGFCtQ85GauWHARzmglpoO3LNsYiY6AxiTGzol8WkwjuaakcvwRoBPnt4XmKKpZYxYcUiy/F7wcgpCo+VfVUqVS/qud56nlEVUFERPPtiUlngIiIJo/BgIiIGAyIiIjBgIiIwGBARERgMCAiIjAYjJyIeJPOw7wQEU9ECiKST3k7qa4/bePOv4jkRKSUxf+FWT8XkjAYDEBEdh1O8JyIVBPWURCRuoioiNSyfHINwx6TmmPyWwB2AayllJeciOwCaMZ8XxSRpv1b5hLWcyAixyJSDi0XnAO7IlK1000RKfRYvmTT1ezPpv0Zu82k/A94jJ0E67T/HwUAI11/dDtprLvPdhPPhUxQVX4cPjAXnbxj2iqAUsL3BQAKwJv0fk3bxx6bwgDpm32OdQlA7Qz5yZl/k8T1a9w2Qn/rZq/5Pc6bznkBIA/goNf+ASi7HKde+Xc9xoMcu+g6R3F+99r+oOfHKD/9zoVRn3vj/vDJwIG9e/dUteW4yG0AlX6JVNU/S76ySFUbqtoYYJGjPt+vniU/jvYAxBWLbNrvXezYnwv259swF5PtmG22B8lkYIBj7Hzshvi7uTi1/ZS2k5ZxnHsjw2DgJiiO6LCPjWX7eN9VLGQv8m0RKbms3K6naYsDSvb3erQYwKbbtcUOVTvPs8UMwXfBo3pQvFC33wdFGrt2mXxc0URonUHxxamii0ieDsJFXnZbxdC6anY9+aT1h4pPdkPrD+9f0/5eiuTJSzhuCzBFd8Vgfr/9s+uqiUgTbsUdNQA+zJ19eD15AO+hf8AK5AD4qhqcO56qbvVKqKptVe0ZDJLyH3OMg2NRE5F6KHnXsQsva7dxbM+jU+sM9seu8yD8dwmfl3baC6WL/ds5nB/R89/p/ypJv3NBEq4D0fw7pJ+sST+azMIH5pG3EJlXtz9zMI/y5cj3NQC7MevrVURwAFPkUQDgAagjVLRgp3Ph9YeWy0XSBXnL23UeACja9QbTZTtdDecllMYL7V/iI79NX45M10PTRdiijn7rjx638P7hpMilFJpXjztudl7TzivZbfXbfldxoF1OE/a9ZI9zOXqcQn+jGmKKidBdJNSEPc9sPo6HOFf75j98jO3fJny8qwnHLmenj+3+VoNt9fi7qZ3nhf4uB5F81kLTwd8hF7f9fudH9Px3+b8awbHseR1IyH/idWOSn4lnYNo/9gTqnKTh+aHfq+ET0M4rx5106B0MohfUThr7Dxtdv2fnH0Tm5yP/VLXIP10JoYtMj3/CUjTf9vtiwjHq7KtdX7TsO7r92PXbZXcjefMiacPT9bjjlrD/Pbdvl40eT5c6g+CCeIyTi2ku9HtiMLDrKEb2q5a03YTzqm/+I8e4aPNdQiTgR49daNlT53V4naFjGr5IJ56XMcvEbT+c937nf+z/1YiOZex1ICb/ideNSX6eBLmKPuoX7KNfXLmtD/OPPsj6/dB0eL3rke+gqr6IfDs6P7RcXN7aCO2LmiKJ8PerMMUuZTvtA9gA0BDTiiNc9LWqph5lD0DVPp5vqmrFFpGUAGzB3AH1XX80ozZvbbvuSrA+PV3XEp7uV46etP0STh/PQdyGKVKshD6J7L70qhNoAqZYQWOKg3rIY8D8q+qePY+qAGr2eK/2OMZhrvkJn2cte565Fpm5cDn/k/6vkrgey37XgbOmHxsGgz7sRRcwJ1cLMM1DYe76lu30AszJE7aM0f3BD2DuVHrNL/eYDwD7Q27LBwDtXVa9aT9dQhfta6HZNZiLuI/uCtSk9Z+iqsu2HLsEoKHulfhxYrdv9yFanuy5rlhVt4Lyd5iyf3/4bOIu7AUa5mbAxcD5t8GmAqBiL1K7OAniQ5EeFenBBTDumPRaxkEa53+g77F0vA4MnX7cWIHspo3utuydP2BwN4yTFiCBTvCIE/kHWED8P+4eer+/cNeuJ3y3fg3AdsKFyOuR17Adu63OP5ljhVsD9rEYMHecdjsbkTvbgdZvn0bqqroVEwiSjls4jcv2GzBPDeHjfCr49bEFc4FyagsfdxG0f7/vAFjrVelpKyKjF5Jh8l8M8mD/Tm10n7dJ50o/4TxXYJ7AAj7MXbJnp2/FrCNp+y7nf+L5YSvAe/2tXI6ly3VgYcD0kzPpcqpZ+MBWloWmPXRXxJZhyl3D7ayb6PFeAswdfh0nFWxBccoxuisQg/LyoNw5XBncKScPzQ+KcEqhbRVt+qACOR/adrgcNchLtPLuwM7vW8kV7FdkXs9l49YfyW9Q4RzUP4Q/QZmx63EL1un12z/7XZB+Fydl0DWcLlMP8lsPbd9Dd9l5kEe1++L1OgcSjqtnlwv/7bvqSXoc29j8R4+x/dTt91WE6oaixy607HHkmEXXGVQ0B5W3nXMrtExQgao4qeDtOs4J2w+fH0nnv8v5UUfyOyJJx9JDwnWgR/4T00/6IzbTlMBG8abaxzuH9DmYC8JMtTOeNva47wKoqH0qsHfDTZzUVxANTUQO0L+OZC6wmMiBPVGq4vjeAMwd1s30cjQ3SjDlzJ2Lvv29jbNV9BLBFhVuMhAYDAaO1LwF6vWr6LKVRDXetY5EtK4meDJoqXsLG6I42zo7bzOnjsVENNXs3dsyQs1T1bElEhG5YzAgIqLZfM/g+eef16WlpUlng4hopjSbzf9U1Yu9vpvJYLC0tIT9/bO+U0JENF9E5OO471iBTEREDAZERMRgQEREGEMwCA0+kTiQg5ixgYtp54eIiE5LNRjYF4SCvtKT0tVg3yq1IwoREdEYpf1k0AZwAQn9uts3dn01w/g1ABxJwjCLREQ0eqk2LQ36/IgMnhK1ju7BT3r1I05ENPd+/r3vAQBe/dnPRr7uaXjPwEP36Ec+evQ/bjuJKwHA5cuXx5AtIqLp4n/wQWrrnpbWRNEBHvxoAlXdVtU1VV27eLHnC3RTbXt7G1tbW9jc3MT6+joajQYajQZWV1exurqKRuOkv6y9vT0sLy9jfX0drVbrVLqtrS0sLy+jUjkpfWu1Wtja2upaR3S9LhqNBvb29rrm+b6PSqXStb2k9EQ0e6YhGPjofhLwMIXjg57F3t4eDg4OUC6XUavVUK1W0W63USgUsLa2hrW1NRQKJ9UkxWIRuVwO6+vryOfzp9KVy2UcHR3h+vXrnWVqtRp2dna61hFdL2Au7Ovr67hw4QKWl5e7gsXm5iZyuRw8z8PqqhmKodVq4cKFC9jePj1Mb6/0RDSbxhkMvPBEqJJ4B91jvOay1q3szs4OlpdPxsXJ5/O4ds0MF+x5HjzPO7VMdH50em1tratLjqOjI7RaLfi+DwBot9s9L9B3795FpVLB8fFx5ykFMHf4nuchl8uhUChgYWEBjUYDuVwOx8fHqFa7WwbHpSei2TSO9wyKMGN9dt4jCEawEhHP9vvfsu8ilJHQ8mhWLSwsoFqtotFodC7W4Qt7MC+s17yw9fV11Ot1AObCvLm5iUKhgLt373bmRZ8KAKBQKHTml8snY4nX6/WugJXL5dBut2ODVVx6IppNqVcgqxkYfS8yz4dpchpMpxIA/t8//RN+d3iYxqo7/ujSJfyXP//zxDS1Wg2VSgWbm5udC+zu7m7nory/v3+qPL7fhTWfz6NWM+N47+7uolarwfd93L59G6VSCc1mE6XS6dc7crmThlqNRgPFonnPz/d9LCycVN14npcYkAZNT0TTbRpaE82FarXaKWrZ3NxEpVJBs2nerysUCqeKYfoFg0KhgHa7Dd/3O3fuxWIRGxsbXfOS1Go13LlzpzN9dHTU9X2/dQyanoimV6aDQb879nFpt9tdd+TVanUkFa65XA43b97sCiSFQgGVSqVTFxCnUqngzp07nQt49M7e9/2uPEcNmp6Ipts0tCbKvKA4JxCU8QPmIhpXZxC92EbT5fP5U4Fmc3MT29vbPesLApVKBbdu3eq6k79+/XqnDgJAp7VTNE+DpCei2ZHpJ4Npsby8jI2NDeRyOSwuLsLzPJTL5c67BkB3he/e3l6npVChUMDR0VHPdOvr66cCRNAsNa7IJggW4XcSyuUyqtUq8vk8KpUKFhcXu5429vb2OnUSe3t7KBaLyOfzsemJaPbM5BjIa2trypHOiGje3P2TPwEAXBvyTWQRaarqWq/vWExEREQMBkRExGBARERgMCAiIjAYEBERGAyIiAgMBkREBAYDIiICg8HYZHWkM9/3sbW1xdHOiGadqs7cZ3V1VWfJ7u6ulsvlznSz2dRaraaqqqVSSUul0qllCoWCVqvVznQ0ned52mw2u77P5/Nd6+i13uPjYy0UCup5nuZyOa3X613pDw4OtF6vd9bVbDYVgHqe17UPwbqKxaIeHx+7HAYiOqOdb35Td775zaGXB7CvMddVPhmMQRZHOgOAjY2Nrp5PiWh2MRiMQRZHOmu1WgBMcNnY2GAxEdGMy3Svpf/nr/8aj95/P9VtLK6s4L/+/d8npsniSGeNRqPTbXWhUMDq6io8z2M31kQzKtPBYJpkbaSzg4ODTnfZQX7q9TqDAdEIPGw2cfzgwan5f3j8GE+cO5fKNjMdDPrdsY9LVkc6i+ZlcXFx4H0gmjdxF/qwzz/6CADwzNJS1/wnzp3Dk+fPp5Iv1hmMQVZHOgvqDQBThxAUORFRvOMHD/DV4WFimmeWlnD5L/4C33jtta7PU5cu4dwzz6SSr0w/GUyLrI50tr6+jq2tLXieh83NTY6BTBQS9wTw1eEhnrp0Cd947bXxZyoBRzojIhrQWYp6fn3vHo4ePBjqDv/R++9jcWUF3713b+BlgeSRzvhkQEQ0oH/90Y9w+O67fStznzx//tRF/zfvvAMAeOHq1YG3u7iyguUbNwZezgWDARHRgD69fx9fHR7i6y+/PPCyL1y9iuUbN/BSj6bfk8RgQEQ0hKcvXx66uGYasTURERHxyYCIZtOvtrdx8OabE9n2F598gqcvX57IttPCYEBEM+ngzTc7rWsG9fjzz/H7L78cettPXbqEr//Znw29/DRiMCCimfT488/xR0O2149r9jmIC9/61tDLTqOxBAMRKQDwVDW2a0ubBqo62GgsRDSXfv/ll/jD48dDLfvM0hIufOtbuDiCbmGyIvVgICI1AFUAOfvCw6mjLyK7ACoAvLg0RERRT5w7N3Vv8s6qVFsT2bt9X1Xb9o7/KHgCiKQ5smlaNo2XZr6IiKhb2k8G6wAOQtNtANEObNoASiJSB9ACUFdVP+V8ERFRSNrBwAMQ7iTft/M6VLUtIhsAdu33G71WJCIlACUAuJyxJl1E8+oszUOz2Lxzksbx0tlCZNoPT4hIDsB1VRWYeoO6iOSjK1HVbVVdU9W1ixcvppZZIhqfoHnoMJ6+fDlzzTsnKe0nAx/dTwIeTLFQWAVADTAXfBFZBbAGU2RERFPuvb/5G3z41ltDLRvc3Q9TCRx0BU2jkfaTwQ5MvUEgFzQdDVUk++iuR1gAwOalRDPiw7fewhcffzzUsme5u3/q0qXMtfWfpFSfDFS1JSItEakCeATzFADbWmhXRK6oakVEaiJShgkMNVVNHgCYiKbK0y++iGsffDDpbNAZpP6egapWeszzAVwITW+mnQ8iIoo3dDGRiLwyyowQEdHkOD0ZiMgVAD/GSdm+ALgC4Gsp5YuIiMbItZioAuBDAK8j0jSUiIhm3yB1BmVV/W1qOSEioolxrTOowr79GxCRH4w+O0RENAmDPBls2SaigKkzUAB/N/osEdEwJjXyF7uFyAanYKCqH4pIRVX/NpgnIjfTyxYRDepXtRr8X/1q7BfmLI76NY+cnwzCgcBO3xl9dohoWL//8ks8dekSVl5/fezb5pvAs8+1aemzAPYAFGCKh/ZU9XqaGSOiwXGwFxqWawXyGwDqMG8NLwDYFpG/TC1XREQ0Vq7B4D1V/VtV/cx+3kb3OAVERDTDXOsMjnvMuw7gpyPMC1EmPGw2cfzgwdi3+4fHj/HEuXNj3y5lg3MwEJEdmLEIPJi6g1pamSKaZccPHkykr/0nzp3Dk+fPj3WblB2uTUvfEpE2zNPAZwD+CsB7aWaMaJY9denSUBW5Z3lX4HeHh1hcWRlqWaLEYCAiK6r6vp1UAP879HUVwP9IKV9EM+vX9+7h0/v38W8/+cnAy/7mnXcAAC9cvTrwsosrK1i+cWPg5YiAhGBgeyp9A8CaiDwHMwxlE+btY8D0WspgQBTx6f37+OKTT4YqJnrh6lUs37iBl0ql/omJRig2GKjqhzBjEUNVPxORdduKCAAgIt8ZQ/6IZtLTly/ju/fuTTobRM6cmpaKyArMsJVhz408N0RENBH96gyehRnQZs1O+/YrD8AtsGkpEVEm9GtNtAxg1/6+gO4XzfZSyREREY1dYjBQ1V8A+GPA1BGE6wyIiCg7XOsMdmBaD4XncXAbIqKMcO2bSFT1jci85VFnhoiIJsM1GLR7zFsbZUaIiGhyXPsm2rFFRXWYlkSbYAUyEVFmOD0Z2Irk12Eqk18G8Lqq3kozY0REND6DDHv5IUxAAHCq3yIiIpphrsNePgfgh+iuJ8gBWEwjU0RENF6uTwa3YDqouwtTV5CHeQmNiIgywLU1EQCUATQAPFLVt2ACAhERZYBrMNgBULIVyW/YF84KrhsRkYKIFPuk8USk3C8dERGN3iCtid6ykxUAzwO45rKsiNRg3lPwRaQZk8YDcAfAtqqyySoR0Zi5ViD/N9gXz6KtivosVwDgq2obQFtEjkSkoKqNSNJdABuq6jvnnIiIRsa1mOh/AfDDM0RkyWG5dQAHoek2TCuk8HqCuodrIrIbV0wkIiUR2ReR/YcPHzpmm4iIXLgGgwqAioi8IiIrdrCbisNyHrq7vfbtvLACTIBo2HXesU8UXVR1W1XXVHXt4sWLjtkmIiIXrk1L9wHUAHw7NM91DORoE1Q/Mr0MYM8WJUFEGjBPFNGiJCIiSolTMFDVz2DHNQg4joHso/tJwMPpTu/8yHQbp4fYJCKiFMUWE4nIPwf1ArZo6JVIkvcc1r8Dc5cfyAWVx6GioB10v7OQBzvBIyIaq6Qngy0ARRFpwVygF0UkKP9fALCBPsVEqtoSkZaIVGHu9itApynprohcsWnqIlKGeUqoBUVGREQ0HrHBQFXfthf/WzAVvB6AVfv1AhzrDFT1VEWzbUJ6ITS9NUimiYhotFzGQL4GACLyp3YadtqlzoCIiGaAc99E4UBgp98efXaIiGgSBumojoiIMorBgIiIGAyIiGiAYCAiSz3eNSAiogxwCgYi8j8BtABUQ/N+kFamiIhovFyfDBZVdQHAdnheCvkhIqIJcA0G/2J/amie80hnREQ03Vx7LRUR+TGAK7Yrietgr6JERJnh2mvpWyLiw/RH9McAfqiqbyUvRUREs8J12Msl+8Yx3zomIsog1zqDDbYeIiLKLtc6gzaAlojchKlEvquqv00vW0RENE7OdQb21zsAICL/ICILqno9tZwREdHYuL509oqIPCsiPxCRf4cZC/luulkjIqJxca0zeBVmFLJlANdUdY2tiYiIssO1zuAAwAVV/SyYISLPst6AiCgbEp8MRGQFAFT1DswLZyvBB6F+ioiIaLbFPhmIyBUAbwBYE5HnYDqqawIQm8RpDGQiIpp+scFAVT8EsGZ//0xE1sNDXXIMZCKi7HBtTbQC4FFk9nMjzw0REU1EYgWyiDwLIAf7hGD7JwIAD8AtAD9NMW9EZ/Kw2cTxgwdj3+4fHj/GE+fOjX27RGfRrzXRMoBd+/sCgKPQd3up5IhoRI4fPMBXh4d46tKlsW73iXPn8OT582PdJtFZJQYDVf0FTC+lEJHvhOsMiGbBU5cu4RuvvTbWbf7bT34y1u0RjYLzGMjRQBA0OyUiotnn2oX1cwB+CFt3ANO89Ao49CURUSa4voF8CyYA3IWpK8jD1CEQEVEGOBcTASjDDHX5yPZLlE8nS0RENG6uTwY7AEqq+ncicldE/gXAuutGRKQAwFPVxBZIIlIDUO+XjsjFr+/dw6f374+9QvfR++9jcWVlrNskOivX8Qx+EXrHoAJgE8B/d1nWXuCrAHIi0lTV1Zh0BQAFAHWX9dLsmFR7/8N3351I09LFlRUs37gx1m0SnVVS30RLMC+Xheet2F93AJQAvJu0cnuB91W1DaAtIkciUlDVRo/keZj+jyhjJtne33vpJXz33r2xbpdoFiU9GWzAPAH4Md+7dFS3DtP9daAN80ZzFxEpwlRMfzthezTD2N6faLolBYMGgO3wGAZhjh3Veeh+a9lH5GnDyqnqnoj0+KqzvRLM0wguX77ssGkiInIV25pIVX8RFwgs147qok1Q/fCEvchv91uJqm7bEdbWLl686LhpIiJy4frS2f+FKRaK+lqfRX10Pwl4MEVFYesANu1TQQ5AXkRqqrrlkjciIjo71/cM9gAsqOrXVPVrAK4D+CuH5XbQ3QQ1F1Qe28plqOqGqq7aVkYNABUGAiKi8XJtWvp6ZHpPRP4ZwJ0+y7VEpCUiVZjxECoAICIegF0RuaKqvp1Xgmla6olIW1XZsojwq+1tHLz55lDLsr0/kTvXYqKlyKxlmJY/falqpcc8H8CFyLxtONQd0Ow5y8tfv3nnHQDAC1evDrws2/sTuXN9A7kBQHEy/vEx7F0+UT+f3r+PLz75ZKj3DF64ehXLN27gpVIphZwRUcA1GFRVNbFIiCjJ05cv8+UvoinmWmdwKhCIyLOq+tvRZ4nSwmEgiSiO65NBMB5yITTruv3QjOAwkEQUx7UC+SZMZ3P7YHcRM43dQhBRL65PBusAVlX1wzQzQ0REk+H60lkNkTeQReSV0WeHiIgmwfXJYB9AU0SOYTqeWwTHQCYiygznpqUwYw3UcBIMimllitLBkb+IKI5zayIAfxluSioi76WQH0rRWV7+Ogu+CUw0/QZ5MqjaPoYCmwBujT5LlCa+/EVEvbgGgyOYdwoKOOmS4gIYDIiIMsH1DeTPbA+jncFuROT76WWLiIjGybVpKXqMeqYjzgsREU1I2iOdUcSk+gcC2EcQEcVLe6Qzigj6B5oE9hFERHFSHemMeptE/0AA+wgionipj3RGRETTjyOdERERRzoDxlup+9E//iP8X/5yIkU27BaCiOLEBgMRuW1/rUUDgYhcAfCnqvrTNDOXpnAA+PyjjwAAzywtpb5d/5e/xJf/8R9j7xICYLcQRBQv6clgVVVf7fWFqn4oImUAMxsM/vVHP8Lhu+92mlo+ef48zj3zTOrb/d3hIS6urrJLCCKaKknBoD22XEzAp/fv46vDQ3z95ZfHul3enRPRNEoKBsd9ls2NMiOTwE7biIiMpJfOLsR9YesMJO57IiKaLUnBYFdE/iE6075z8HMAP04rU0RENF6xxUSq+raILIvII5j3DI5gioYKAPZmuSURERF1S3zPQFW3RaQNM8TlMkyl8quq+vY4MkdEROPR96UzVW3APBkQEVFGOY9nQERE2TWWYCAiBREpjmNbREQ0uNSDgYjUYOoafBFp9vjeE5G6iByLyIGIFNLOExERdUs1GNgLu6+qbVv3cNTjYn8NpiO8CwBqAOpp5omIiE5L+8lgHcBBaLqN028uN2yggKpupZwfIiLqIe1g4MG8nxDw7bwOVe30gWSfGvZ6rUhESiKyLyL7Dx8+HHlGiYjm2TgqkBci035C2k0AN3t9oarbqrqmqmsXL14cVd6IiAjpBwMf3U8CHmJ6QxWRKoCbquqnnCciIopIOxjswNQbBHJB/UC4ItkGgtsMBEREk5FqMFDVFoCWiFTtYDgVwDQnhekIz7NNT8sAjkVE7aeaZr6IiKib6xjIQ1PVSo95Pk66yN60HyIimhB2R0FERAwGRETEYEBERGAwICIiMBgQEREYDIiICAwGREQEBgMiIgKDARERgcGAiIjAYEBERGAwICIiMBgQEREYDIiICAwGREQEBgMiIgKDARERgcGAiIjAYEBERGAwICIiMBgQEREYDIiICAwGREQEBgMiIgKDARERgcGAiIjAYEBERGAwICIiMBgQERHGFAxEpCAixbOmISKidKQeDESkBqANwBeR5rBpiIgoPakGAxEpAPBVta2qDQBHdt5AaYiIKF1Pprz+dQAHoek2gNwQaUbi59/7HvwPPgAAfPHxx3j6xRfT2AwR0cxJu5jIA3AUmvbtvEHTQERKIrIvIvsPHz48c8aefvFFXPn+98+8HiKiLEj7yQAAFiLT/jBpVHUbwDYArK2t6TAZefVnPxtmMSKizEv7ycBH912+B1MMNGgaIiJKUdrBYAemTiCQs5XECFUSx6YhIqLxSLWYSFVbItISkSqARwAqACAiHoBdEbkSl4aIiMYn9ToDVT11cVdVH8CFpDRERDQ+7I6CiIgYDIiIiMGAiIjAYEBERABEdaj3tyZKRB4C+HgEq3oewH+OYD2zZN72mfubbfO2v8DZ9vlFVb3Y64uZDAajIiL7qro26XyM07ztM/c32+Ztf4H09pnFRERExGBAREQMBtuTzsAEzNs+c3+zbd72F0hpn+e6zoCIiIx5fzIgIiIwGBARTRUR8USkajvvjH5XEJGi6/xBzG0wGMXBm2VZ3n8RqUX3Lav7ay8c5Tna30KvMdKzsr8ikgdwDKDU47sazFgvvog0+80f1FwGg1EdvGlmLxJ1ETkWkYPwP1CW99/uZyEyL5P7a7uCvwNgW1X3QvOzur+7MPt1lMbFcEq0YXp07urJ2Z7Xvqq27XgvR6HAeGr+MBueu2AwyoM35a4BqKrqBQA1AHVgLvY/D6AVTGR8f3cB3LRdwgPI7v7afTiy+9WC2S8va/urqn747xmyDuAgNN0GkEuYP7C5CwYY4cGbco1gxDhV3QrNz+z+22KC4A7Ztz8zub+2OAEAronIbqiIJJP7C7MfJREpikgOQN1eNLO6v1EegKPQtG/nxc0fWOqD20whDyM6eNNMVTvjSNs7peAi6SG7+59T1T0RCc/zkM39LcBc9Br20xQRHxndX1Vti8gGzNOQD2DDfuUhg/sbYyEy7feZP5B5DAbAiA7eDNkEcDM0nbn9F5ES4l/Gydz+AlgGsBcEfRFp4GQs8cztr30auK6qYv/WdRFZtV9nbn978NEd5DyYp6DlmPkDm8diIh8jOnizwDZPC5cr+8jm/q8DeNtWIBYA1ESkjOzurx+ZbsOMIe4jm/tbgan7gqpuwwT+NWR3f4Hu/drBSbAHzFNwI2H+wOYxGIzs4E07GwhuRyqkMrn/qrqhqququgpTbFKxdSWZ3F+Y/cqHpvMwRYFZ3V8f3XUBCzB/58ztr63/2QTQaS5rK81b9v2DMmxro7j5Q213HrujCL3M8QimorWVlH4W2eZ20bbKW6payfL+2yKEKoB9mIDQyur+2n9+wFwoj4LmpRne3xpMZbEPIGg9lNn9Hbe5DAZERNRtHouJiIgogsGAiIgYDIiIiMGAiIjAYEBERGAwoDlj+7Y5sD1g9vq+ant6LdleIQ9EpGm7ia7Z3/Oh9AU7v2p/1nv1Qx/ZRt7+LNlPMfRd0fZG2pWWKG3z2h0FzSnbd9EmgKKIeD16iCzBtNnfBjrdPPhBZ3/2Qn8HwKp9p2FVVTfDK0gKBratfMUGAE9Vt2yAadi85CJ5aotILboNolHjkwHNozbMi0u3wjPtxX0/ktaPTL9n0xZgugjvdZGu9dqoXX/TXuy/jZNuE44ALER6XQVgujQGcJCFgVtoujEY0DzyAdzG6Te0NxBzIQ9Zh+kCYQOmO4RTwj3GRlSCJw67juu2AzbfLpOLWXYbkcBFNGoMBjSXQsU+JfszDzMAkN8rva0bKOOkyCiH7q6TXXR617RdJtyEKSrasE8accHFRzb76KcpwmBA82wbJx17bUYGAeqiqg1V3VLVIP1Ag6jYYNN1129HtQr60cnbfpTKtlK5Fq5IhhnBiwGBUsNgQPPsNoCcfTrwB1x2F6ZXyUEu0F6vmTZQtOzFf90WJTXRY1B0orQwGNDcssUvWzD1BLftbC+SzOsxD7bHzD2YQVa6AkJMgEh6kijY9YUHafHQHaDi6hOIRoJNS2mu2FY5RRE5sHfgt2HK7X17Z34dJ08LbZiBciAiwQW7w5b1F2EG0gFOioHq6FEkJCKtaHNWGzjaNk3bpikCWA5aKtk0M91HP00/dmFNNCb2Ip8L103EvOsQXa4KMwA8AwKlhsVERGNiB59ZDlcMOwQCz6ZjIKBU8cmAaMxEJO86GtcgaYnOgsGAiIhYTERERAwGREQEBgMiIgKDARERgcGAiIjAYEBERAD+Pxt62G11nj7qAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "set_texfig()\n",
    "plt.step(stats2010['MPC'], stats2010['cdf_weight'], where='post', color='#9A0000', alpha = 0.5, label='SHIW 2010')\n",
    "plt.step(stats2016['MPC'], stats2016['cdf_weight'], where='post', color='#9A0000', label='SHIW 2016')\n",
    "plt.xlabel('MPC (\\%)')\n",
    "plt.ylabel('Cumulative fraction')\n",
    "plt.legend(framealpha=0)\n",
    "plt.title('(a) Income-weighted MPC distribution, data')\n",
    "plt.savefig('figures/figC1_a.pdf', format='pdf', transparent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6738d03",
   "metadata": {},
   "source": [
    "# Figure 1 asset response lower bound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a87cfe48",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = calibration.r  # 0.05\n",
    "T = 6\n",
    "\n",
    "y_weight = (s2016['y'] * s2016['pesofit']).values / (s2016['y'] * s2016['pesofit']).sum()\n",
    "mps = 1 - s2016['riscons2'].values/100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07ac7b4d",
   "metadata": {},
   "source": [
    "Implement equation (A.48), which gives an upper bound on asset holdings in response to a date-0 income shock, in terms of the income-weighted MPS distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "96708807",
   "metadata": {},
   "outputs": [],
   "source": [
    "dAubar = np.array([(1+r)**t * y_weight @ mps**(t+1) for t in range(T)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee516f51",
   "metadata": {},
   "source": [
    "By corollary 2, if we calculate $\\underline{dC_0} \\equiv 1-\\overline{dA_0}$ and $\\underline{dC_t} \\equiv (1+r)\\overline{dA_{t-1}} - \\overline{dA_t}$, then the cumulative PDV of $\\underline{dC_t}$ is a lower bound for the cumulative consumption response up to any point.\n",
    "\n",
    "Calculating, we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "625050e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "dClbar = - dAubar\n",
    "dClbar[0] += 1\n",
    "dClbar[1:] += (1+r)*dAubar[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0ff13c7c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD6CAYAAABTcqc2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnGElEQVR4nO3de3xU9Z3/8dd3kpAbkCEEAggIAQXkJhC8oXghsbW1tlXQ3rS2W4Nrt9a6bVK2/VXZ3S4btleta5Nf++tFd1sMq9ubNwIoVK0CUVEEFIKKKCGQDBIuuc3398eZCZP7hZmcubyfj8c8mHPmzMwnKO/zzff7Pd9jrLWIiEhs87hdgIiInDmFuYhIHFCYi4jEAYW5iEgcUJiLiMSBZDe+NCcnx06aNMmNrxYRiVnbtm07bK0d1dVrroT5pEmT2Lp1qxtfLSISs4wx73T3mrpZRETigMJcRCQOKMxFROKAwlxEJA64MgAqIonJ7/fz3nvvcfz4cbdLiVqZmZmMHz8ej6d/bW2FuYgMmsOHD2OMYdq0af0Oq0Tg9/s5cOAAhw8fZvTo0f16r/42RWTQ+Hw+cnNzFeTd8Hg85ObmcvTo0f6/NwL1RMyRhkYefGYvza1+t0sRkQFobW0lJSXF7TKiWkpKCi0tLf1+X0yF+V9e+4DSJ3dx3c+e4/UD/T9ziYj7jDFulxDVBvr3E1NhfsvFkyi/eQGHGxr55APP8YOndtPY0up2WSIiroupMAe4euYYKr9xOZ+edxY/27iHj9/3V15+t97tskREXBVzYQ6QlZHCD5bN5ddfWsjxxhZuePB5vv+XNzjZpFa6iESXqqoqfD5fxL8nJsM86Ippo3n6G4v5zAUT+b+b93HNTzfx0r46t8sSEWmzbNkyqqurI/49MT/PfFhaCv/26dlcO3ssJY9u58ayF/jixWdT/NHpZKbG/I8nEtdW/mkHb7z/YUS/47xxw7nnEzN7Pa6qqoo1a9YAsGLFCrxeLz6fj/LycrxeL+vWraO0tJTs7OxO+/Ly8li7di0AW7ZsYeTIkRQXF7d9dnCV2Ly8vC7fGw4x3TIPdcnUHJ78+mJuvWQSv/3bO3zkJ5v461uH3S5LRGJAdXU1lZWVrFixAp/PR0lJCUBb8BYVFXUK8tB9lZWVrFu3jqVLl1JaWsqqVavaPrOu7nRvQVfvDZe4arpmpiZz73Uz+ficsZSs3c4Xfvkin71gAis+NoPhaZrbKhJt+tJiHgx5eXkUFxdTVVXFlClT2lro8+fPp7CwkHXr1nHTTTexdOnSLvdVVFRQXV1NeXl52+dVV1dTUFBAdnY2BQUF5OXlUVdX1+m94RI3LfNQCydl8/jXL2P55Xms2bKfq3+0iY27DrldlohEqerqapYtW0ZeXl67gC0oKGDv3r0sXLiQkpISVq9e3eU+gMLCQoqKiigqKmLbtm0UFBR0+p7u3hsOcRnmAGkpSay4ZgaP3bGI4enJfOnXW7j7kVfwnWhyuzQRiTJlZWUsXLgQr9dLVVVV2/7y8vK2VntpaWm3+woLC9ta810Jzmbp6r3hElfdLF2ZO8HLn752KQ9s2MN/PrOXTW8e5l8/NYuPzhrjdmkiEiVuuukmli1bxt69e1mwYAHV1dVtUwqXL1/OlClTACguLmb16tWd9oEz8LlgwQIKCgraDYAuXbqU2267ra0/vqv3hoOx1obtw/oqPz/funEP0B3vH6V47XZ2vP8hH58zlpXXzSRnaOqg1yGSqHbu3MmMGTPcLiPqdff3ZIzZZq3N7+o9cdvN0pWZ47L4368u4ptXn8u6HTVc/eNN/PHV93HjhCYiEk4JFeYAKUke/uGqc/jznZcyITuDO3/3MkUPbePQh6fcLk1EZMASLsyDzs0dxqN/fwnf+dgMNr1ZS8GPnqVi63610kUkJiVsmAMkeQy3Lc7jybsWM33McL61dju3/moLB3wn3S5NRKRfEjrMgybnZPL7ootYed1Mtrxdx0d+vIn/evEd/H610kUkNvQpzI0xBcaYbi9VMsZUdHh4w1bhIPF4DF+8ZBJP3bWYuROy+M5jr/P5X7zIu0dOuF2aiEiveg1zY0wZUA34jDHbejh0TfBhrfWFp7zBNyE7g4f/7kJWXT+b1w4c5SM/2cSvntunVrqIRLUew9wYUwD4rLXV1tpKoC6wr6Mt1tq1wUdEKh1Exhg+e8FEnv7GYi7My2bln97gxrIX2Fvb4HZpIhIBa9euZcqUKSxbtqzL10tKShgxYgTl5eVUVlYyZcoUFixY0HYB0YIFC9pdOVpZWcny5cspKSlh+fLlFBYWti3eFTHW2m4fQClQFLJdFrod2OcFKoB1gcfSbj6rCNgKbJ04caKNFX6/3/7Ptv12zr1P2XO/87j9+TN7bHNLq9tlicSkN954w+0SulVQUGABW19f3+k1r9dr8/Ly2raLiopscXFx23ZxcbGdP3++tdbasrIyW1RU1OkzQo/vTXd/T8BW201e99bN4gVC7/bgC+zraIu1tjAQ9hXGmE4L9Fpry621+dba/FGjRvXhNBMdjDFcP3886+5ezBXTRrHqiV3c8ODz7D54zO3SRCSM8vLy8Hq9rFq1qt3+8vJy8vPbX3Tp9XrbbS9cuBBwWuQlJSWUlZV1+vzly5eHt+AO+jIA2nHBXV/ohrXWZ61dHXi+FqgEuuqKiWmjh6Xx8y8s4Gefm8f++pNce/9m7lv/Fs2tfrdLE5Ew8Hq9rFixom0Z26CKiopegzi4pG1FRUWXqyUCYbsJRXd6W2jLR/uWuBdnMLQn1X04JiYZY7h2zjguzhvJyj+9wY/WvckTrx/kP5bOYdZZWW6XJxJ7nvg2HHwtst8xZjZc8+99OrS4uJhVq1ZRXl5OUVERVVVVFBYWdmqJB1VWVlJVVYXX66W4uJjCwsKIh3Z3emuZrwEKQ7bzrDMQGhwcxRjTadmv4DHxauTQVO777DzKb17A4YZGPvnAc/zgqd00tuiG0iKxLngXIHCWxu1pZcOCgoJ2y9kGb0rhhh5b5tbaKmNMlTGmFDgClAAE5pFXGGMmA1XGmApgC87Ml8h2DEWRq2eO4cLJI/mXv7zBzzbu4akdB1m9dA7zJo5wuzSR2NDHFvNgWrFiBatXr267xVt/LFu2jMLCQqqrqwe9hd7reubW2k7zaawzjzyYWJWBR0LKykjhB8vmcu2csax49DVuePB5vnJZHncXnktaSpLb5YlIPwW7TJYvX059fT1w+uYSQR23gwoKCli6dGnbreFCAz3SAR/3N6cYLFdMG83T31jMqid2Ub6pmnVv1FB6wxwumBy+G7aKSGSsXbu2ba55UVFR240kvF4vPp+PNWvWtN3jMy8vj8pKp/1aWVnZacCzoqKCtWvXtg2aBgM80v3pCXVzisHy/J7DlDy6nf11J/nixWdT/NHpZKbqvCmim1P0jW5OESUumZrDk19fzK2XTOK3f3uHj/xkE8/tOex2WSISxxTmEZKZmsy9183kkeUXMyTJw+d/8SIrHt3Oh6ea3S5NROKQwjzCFk7K5vGvX8byy/NYs2U/H/nxJjbuOuR2WSISZxTmgyAtJYkV18zgsTsWMSwtmS/9egt3P/IKvhNNbpcmInFCYT6I5k7w8qevXcqdV03lj6+8T8GPNvHk6wfdLktkULkx6SKWDPTvR2E+yFKTk7j76mn84R8WkTs8ldsf3sZX/7uKIw2NbpcmEnFJSUk0N2vcqCfNzc0kJ/d/9pvC3CUzx2Xxv19dxDevPpd1O2oo/PEm/vjq+2q1SFzzer3U1NTg92uBuq74/X5qamrIyur/Wk+aZx4F3qw5xrfWbufV/T4Kz8vl+5+axejhaW6XJRJ2fr+f9957j+PHj7tdStTKzMxk/PjxeDyd29o9zTNXmEeJVr/l//11Hz94ejepyR6+94mZ3DD/LIwxbpcmIlFCFw3FgCSP4bbFeTx512KmjxnONyte5dZfbeF930m3SxORGBBbYd50AnY85nYVETU5J5PfF13EyutmsuXtOq7+8SZ++dd91B7TAKmIdC+2ulleeACe+if40pNw9sXhLyzK7K87wbcf3c5ze44AztTGgumjWTIjlxljh6kLRiTBxE+fedMJuH8+eCfCl5+CBAgzay07PzjGhl01VO48xKvv+bAWxmWlcdUMJ9gvzhup5XZFEkD8hDnA1l/Bn++Cz/wOpn8srHXFgtpjjWzcdYj1u2rY/NZhTjS1kp6SxKXn5FAwYzRXTh/N6GGaCSMSj+IrzFtb4D8vBE8y/P3z4EncFump5lb+Vn2E9TsPsX5nDe8fPQXA3PFZLJmRy5IZozlv7HB1x4jEifgKc4A3/gCP3AKffADmfSF8hcUway27Dh5j/c4a1u86xCv7ne6YsVlpXDV9NAUzcrl4irpjRGJZ/IW5tfCLJXDsIHxtG6Skh6+4OFF7rJGNuw+xYechNr1V29Yds2iq0x1z1fTRujBJJMbEX5gD7NsMv7kWCv8FFt0ZnsLiVGNLK3+rrmPDTmcQ9UBg7vqc8Vksme50x8wcp+4YkWgXn2EO8PBSeG8LfP1VSPee+eclAGstu2uOtfWzvxzojhkz3JkdUzBjNJdMyVF3jEgUit8wP/ga/PwyuPQuKLj3zD8vAR1uCMyO2XmIzW/VcryplbQUD5dOzWHJjFyumj6aXHXHiESF+A1zgEeLnAHRO1+G4ePC85kJqrGllRer61jfoTtm9llZLJnhDKKqO0bEPfEd5vVvw/35cP7n4Lr7wvOZgrWWN2saqNxZw4Zdh6h6tx5rIXd4KldNz23rjkkfou4YkcES32EO8MS34aUyuONFGHVu+D5X2hxpaGTj7lo27Kph05uHaWhsIS3Fw6IpOW1z2tUdIxJZ8R/mxw/DT8+HKVfATQ+H73OlS00tfl7c51ysVLmzhvfqne6YWWcNZ8n03LbuGI9H3TEi4XTGYW6MKQC81tq1vRxXBqzr7biIrGf+7GrY+H34u0qYsDC8ny3dstby1iGnO2b9ztPdMaOHpbJkxmiWTM9l0VR1x4iEwxmFeSCgS4E8oNRau6Cb4wqAMqDElTBvbID75kHOOXDrXxJiEa5oVHe8iY27DrFh1yGefbOWhsYWUpM9LJqaw5LAxUpjs3SRl8hA9BTmPd41NBDQPmttNVBtjKkzxhRYayu7OHw+UHXm5Q5Q6lC4vBge/ya8tQ7Ovdq1UhJZduYQblgwnhsWjKepxc9L++qcVvsuZyAVYOa44U4/+/TRzD4rS90xImHQY8vcGFMK7LXWlge2y4Btwe2Q45biBHkpUNZN2LeJ2G3jWpvhZwshJQNu35zQi3BFG2stew41ULnzEBt21bDtnXr8FkYNS2XJdKfFfuk5OWQM6f9dyUUSxYBb5oAXqAvZ9gX2dZRnrV3b0/xjY0wRUAQwceLEXr52gJJSYMn/gbVfhtcqYO5nIvM90m/GGM7JHcY5ucP4+yumUHe8iWd2Oxcr/WX7B/x+y35Skz1cMmWks0b7lJGMH5FOarJOyCJ90ZdmUHaHbV/oRiCky+lFoDVfDk7LvI/19d95n4axP4UN34eZn4bk1Ih9lQxcduYQrp8/nuvnO90xW96uaxtE3bj7dcAZ9sgdlsaE7HQmjMhgQnbgMSKdCdkZ5A5PI0ldNCJA72Huo31L3AtUdzimEFgeaJXnAfONMWXW2tXhKbGfPB4oWAkPfQq2/BIuvsOVMqTvhgQGSBdNzeF7157H3toGXt1/lP31J9hfd5L9dSd4ofoIj71ygNBewZQkw1leJ9jHj8hgQnY6E7Mz2oJ/REaKrlaVhNFbn/l8nBkshYHtdSHPOw2EGmMqgDWuzGbp6LefdNZuufMVSBse2e+SQdHY0sr7vlPsrztxOujrTzjbdSeoP9Hc7vjMIUntgn7CiAwn7LOdbfXPS6wZcJ+5tbbKGFMVGAg9ApQEPtALVBhjJltrfYF9RUAB4DXGVFtr3ZvZAs7CW+VXwPP3w1XfcbUUCY/U5CQm52QyOSezy9cbGlvagn1/vdOify8Q9s/tOczJ5tZ2x4/MHML4kG6b02GfzjhvOilJnsH4sUTCIj6uAO1OxZfgzSed1vmw3Mh/n0Qtay1Hjjd1EfRO6/5A/Ula/Kf/LXgMjM1KZ3xI0E/Idp5PzM5g1NBUTamUQXcms1li21XfhZ1/hE2r4eM/dLsacZExhpyhqeQMTWXexBGdXm/1Wz44erIt3N8LCf3Nb9VS82Fju+OHJHucoA/pwglt3WdlpAzWjyYCxHuYj5wCC26Fbb+Gi+5wtkW6kOQxjB/h9K9fzMhOr59qbuW9+s5Bv7/+BK/s93H0ZPv++mFpyZ2DPvB8/IgMLW8gYRff3SwAx2qcy/zP/Qgs+9XgfKcknKMnmzt13YR26TS2+NsdP2pYaru++gnZ6eQMTSU9JYm0IUmkpwQeQ5JICzxPSTKanZPgErebBZy+8ou/6nS1LLoTxs1zuyKJQ1npKWSdlcWss7I6vWatpbahMTA4e7LdbJxt79Tz5+0f0OrvvVGV5DFkdAh757mnU/B32g68Jy3kubPfQ1pKEhlDkklPSSI12aOxgBgV/y1zgFMfwn3nw5jZcMsfBu97RfqgudXPwaOnqD/RxMmmVk42t3Kq2fnzZJP/9HZTKydCXw8877jddmxzKwP5550WPDn08FtC55OFp5uTRfvt4Of1NFPIWovfOuMYfmtp9Vta/Ba/39JqT//Z6j/9cI47/Z6WdvtPv6ftc0LfYy2tfj+tftp9tt9aWlpPf0bbd/txjg+839+hltPHnX7eElLLp84/i2tmjx3Q/yuJ3TIHZ5754m/Bk9+GvRtgylVuVyTSJiXJ03Z1azhZa2lq9XMqcEI42U3gnwo5KZxsCj2RtD+xHG9s4XBDU6cTSVOHLqS+SPYY0lOSwBAS0tDi99OHX1KiQpLHkGQMHg8kGeNsBx4e0/7PZI/BEzi+4/hKuCRGmAPkfxn+9p9QeS9MvsK5UlQkjhljSE1OIjU5iSwiN7um1W87nQDa/dbQzcniRJPzm0NoCDrhaALhCJ5gEJrOQZnUdlz79yR5PM57TfvPTG53XIfPCQ3mkGOSunhP8HOiTeKEeXIqXPldeKwIdjwKs5e6XZFIXEjyGDJTk8lMTZw4iUaJ1TydvQxyZ8GGf4WWJrerEREJm8QKc4/Hucy/fh9U/cbtakREwiaxwhxgagGcfSk8W+rcak5EJA4kXpgbA4Ur4XgtvPCA29WIiIRF4oU5wPh8mPEJeP4+OH7Y7WpERM5YYoY5wJJ7oPkkbPoPtysRETljiRvmOefAvC84dyOqf9vtakREzkjihjnAFd8GT7Jzv1ARkRiW2GE+fBxcdDu8VgEfbHe7GhGRAUvsMAdYdBekZcH6lW5XIiIyYArzdC9c9o+wpxL2bXK7GhGRAVGYA1xQBMPHO4twubAksIjImVKYA6SkwZUr4MA2556hIiIxRmEeNPezMGoGrP9naG1xuxoRkX5RmAd5kmDJ9+DIHnj5IberERHpF4V5qGnXwISL4Jl/h6YTblcjItJnCvNQwUW4Gg7Ciw+6XY2ISJ8pzDuaeBFM+xj89Sdwos7takRE+kRh3pUl34OmBtj8Q7crERHpkz6FuTGmwBjT7U0zjTFFxph1xphSY0xe+MpzyegZMPdz8FI5+Pa7XY2ISK96DXNjTBlQDfiMMdu6eL0AmAIsA7xARZhrdMeVKwADz6xyuxIRkV71GOaBoPZZa6uttZVAXWBfO9baEmutDyjBCfTYlzUeLiyCV/4bat5wuxoRkR711jIvBPaGbFcD7bpRAiEflA+Uhae0KHDp3ZA63LmQSEQkivUW5l4gdEqHjy5a3sYYrzGmGCgFKju+HjimyBiz1Riztba2dkDFDrqMbLj0LnjzCXjnBberERHpVl8GQLM7bPs6HhDoYinH6WbZZoyZ38Ux5dbafGtt/qhRowZQqksuvB2GjYXKe7QIl4hErd7C3Ef7lrgXp6ulE2utL9DlUonT3RIfhmTA5SWw/0XY/bjb1YiIdKm3MF+D028elBfsIw8OhAa6V0LlAY+ErcJoMO9mGDlVi3CJSNTqMcyttVVAVWD+eDFONwrGGC9QEfizOjDHvNgYUwQsCHS7xI+kZOdCotpd8Orv3K5GRKQTY13oB87Pz7dbt24d9O89I9bCL5bAsYPwtW2Qku52RSKSYIwx26y1XXZj63L+vjIGClbChwecK0NFRKKIwrw/Jl8GUwth84/gZL3b1YiItFGY91fBPXDqqLOqoohIlFCY99eY2TDnRnjx5/Dh+25XIyICKMwH5sp/An+rFuESkaihMB+IEZNg4Vfg5Yeh9k23qxERUZgP2OJvQkomrF/pdiUiIgrzAcvMgUV3wq4/w/4tblcjIglOYX4mLroDMkdrES4RcZ3C/EykDoXLi+Gd5+CtdW5XIyIJTGF+phbcCiMmQ+W9zgwXEREXKMzPVFIKLPk/cGgHvBYftz8VkdijMA+H8z4NY+fChu9DS6Pb1YhIAlKYh4PH4yzCdfRd2PJLt6sRkQSkMA+XKVdC3hWw6T+ctVtERAaRwjycCu6Fk3Xw/P1uVyIiCUZhHk7j5sHM6+GFB+BYjdvViEgCUZiH21XfhdYmeLbU7UpEJIEozMNt5BRn7nnVb+DIXrerEZEEoTCPhMXFkJQKG/7F7UpEJEEozCNhWC5c/FXY8RgcqHK7GhFJAArzSLnka5Ax0rnMX0QkwhTmkZI2HBZ/C/Y9C3s3uF2NiMQ5hXkk5X8ZvBNh3T3g97tdjYjEMYV5JCWnwpXfhYPbYcejblcjInFMYR5ps5dB7ixnZktLk9vViEicUphHmsfjXOZf/7Yz91xEJAL6FObGmAJjzNJIFxO3phbA2Zc6V4U2NrhdjYjEoV7D3BhTBlQDPmPMti5e9xpj1hlj6o0xe40xBZEoNKYZA4Ur4Xits26LiEiY9RjmgWD2WWurrbWVQF0XYX0jUGqtHQGUAboZZlfG58OMT8Dz90FDrdvViEic6a1lXgiELjBSDeR1OKYyEPRYa1eHsbb4s+QeaD4Jm3/gdiUiEmd6C3MvUBey7Qvsa2OtrQ4+D7Ta13b1QcaYImPMVmPM1traBG2Z5pwD877g3I2obp/b1YhIHOnLAGh2h21fD8cuB27r6gVrbbm1Nt9amz9q1Kg+lheHrvg2eJJh47+5XYmIxJHewtxH+5a4F6erpRNjTClwm7XWF4a64tfwcXDR7fDaI/DBdrerEZE40VuYr8HpNw/KC/aPhw6EBoJ8lYK8jxbdBWleWL/S7UpEJE70GObW2iqgyhhTaowpBkrAmY4IVASmJZYBxUC9McYGHrrNTk/SvXDZP8KeSti3ye1qRCQOGGvtoH9pfn6+3bp166B/b1RpPgX3L4Cho+G2Dc5cdBGRHhhjtllr87t6TZfzuyUlDa5cAe9XwRt/cLsaEYlxCnM3zf0sjJrhLMLV2ux2NSISwxTmbvIkwZLvwZE98PJDblcjIjFMYe62adfAhIvgmVJoOu52NSISoxTmbgsuwtVwEP72oNvViEiMUphHg4kXwbSPwXM/hRN1vR8vItKBwjxaLPkeNDXA5h+6XYmIxCCFebQYPQPmfg5eKgffu25XIyIxRmEeTa5cARjYuMrtSkQkxijMo0nWeLiwCF79HdTscLsaEYkhCvNoc+ndkDoc1v+z25WISAxRmEebjGy49C5480l453m3qxGRGKEwj0YX3g7DxsK6e8CFhdBEJPYozKPRkAy4vATeewl2P+52NSISAxTm0WrezTByKlSuhNYWt6sRkSinMI9WScnOhUSHdzuzW0REepDsdgHSgxnXwVkL4OnvwtH9cP7nYcTZblclIlFILfNoZgx86kEn0J9dDT+dCw99GnY8Bi2NblcnIlFELfNoN2oa3Pyoc4n/y/8FLz8MFbdCxkjn5hbzbobR092uUkRcpnuAxhp/K1RvhKrfwq7Hwd8MEy6E+bfAeZ+C1KFuVygiEdLTPUAV5rGsoRa2/94J9sNvwpBhMPsGJ9jHzddNokXijMI83lkL+190Qv31R6HlJOTOckJ99jLnqlIRiXkK80Ry6ii8/j9OsL//MiSlwoxPOME+6TLwaMxbJFYpzBPVB9udG0VvX+OE/IhJzoDp+Z+H4WPdrk5E+klhnuiaT8LOPzmt9bc3g0mCc652WuvnXO1coCQiUa+nMO/Tv2JjTAHgtdau7eZ1L7ACwFpbMsA6JVJS0mHOjc7jyF6ntf7Kf8ObT8DQMXD+52DeF2DkFLcrFZEB6rUD1RhTBlQDPmPMti5enw/UA0XhL0/CbuQUKLgXvrEDPvM7GDcPnvsJ3D8ffn0tbK+A5lNuVyki/dRjyzzQIvdZa6uBamNMnTGmwFpbGXJYNTACuBFQ0y5WJKXA9I85jw/fd1rqLz8Ej34F0rww5yanG2bMLLcrFZE+6K2bpRDYG7JdDeSFHmCt9QEYzWmOXcPHweJvOnc5enuz07e+7VfwUpkzX33+LTDrBkgb7nalItKN3sLcC9SFbPsC+yQeeTyQd7nzOFEH2x+Bqt/An++Cp/4JZl7vBPuEC3RBkkiU6csAaMcrTnwD+SJjTBGBfvWJEycO5CNkMGVkw0W3w4XL4UCVE+qv/w+88jDkTHNCfe5nIDPH7UpFhN4HQH20b4l7cbpa+s1aW26tzbfW5o8aNWogHyFuMAbGL4Dr7oN/3A3X/QzSvfD0d+CH0+GRW2BPpbNmjIi4prcwX4PTbx6UFxz8DAyOduQNU10SjVKHwvyb4e+ehjtedFrt+zbDwzc4y/M+8+/g2+92lSIJqdeLhowxpYGnR4BKa21VYF75PmCytdZnjFmKM8/cC5R0Nx89SBcNxZGWRuc+pVW/hb0bnX1TlzjdMOdeA8lD3K1PJI7oClAZHPXvwCuBNdc/PAAZOXD+Z2HeLTDqXLerE4l5CnMZXP5W2LvBGTTd/QT4W2DixYE11z8JQzLdrlAkJinMxT0Nh5wbUlf9Fo7sgdThMHupE+xjz9cUR5F+UJiL+6yFd19wQn3H/zprro+ZDfO/6IR7+gi3KxSJegpziS4nffD6Wqh6CD54BZLTYMZ1sOCLcPYitdZFuqEwl+j1watOqG9/BBqPQnaes+b6rBvAO1HBLhJCYS7Rr+nE6TXX3/mrsy/N63TF5M5y/hwzC0ZNh+RUV0sVccsZr2cuEnFDMmDuTc7j8B7Y9wwcfB0OvubMimk+4RznSXaWExgzq33Qa1kBSXAKc4k+OVOdR5C/FeqqnWCvCQT8vs3O7fCCho1t34IfM8fpsvEkDX79Ii5QmEv08yRBzjnOY9b1p/cfPwI1r51uwde8DtUbnXntACkZMHpGSAt+DuSeB6nD3Pk5RCJIYS6xK3Mk5F3hPIJaGqF29+kW/MHXnKmQ2359+pgRkwMt+Nmngz5rvAZbJaYpzCW+JKfC2DnOI8haZ3mBthZ8IOR3/vH0MRpslRinMJf4Z4zT8s4aD9M+enp/YwMcegMObtdgq8Q8hbkkrtShzl2TJlxwep+/Fer2OQGvwVaJIQpzkVCepNOzaTTYKjFEYS7SF30dbH3jDxpsFVcozEUGqs+Dra87V7cSuNpag60SAQpzkXDq62BrzeudB1uzJsCwMTA0N/DnaBg6BoblOn8OzYWMkeDp7W6PkogU5iKDoS+DrfXvQEMN1Oxwbu7R+GHnz/EkQ+ZoJ+i7DP7g81y19BOMwlzELd0NtgY1nXDCvaEGjh10bvTRcBCOBfZ9eAAOVMHxWtq6cEKljwi06DsGf27756nD1IcfBxTmItFqSAZkT3YePWltgROHA4HfMfgD+959wTkJtDZ2fn9yekhXTg/BnzFS0y+jmMJcJNYlJTthO2xMz8dZC6d8p1v2bcEf8rx2F+x7Fk4d7fx+kwSZo3oP/qG5kJIWkR9VuqcwF0kUxjhdL+kjYPT0no9tPhkI+I7BH2j1H/vAuUvU8Vqw/s7vT8vqMHgbDP7gvsAjLUtdPGGiMBeRzlLSYcQk59ETfyscP9y+L7/teSD497/o7G851fn9xgNDhjn99qnDnIHi4PN2+4OvDQ+8FnJc6nDntQQf8FWYi8jAeZKclvawXBjbw3HWOl037QZxDzr7Go+1f5z6EI4egKaG0/u6GuDtKGlINyeCDsHf7kQQcoIIPUnE4NiAwlxEIs8YSPc6j1Hn9u+9fj80H3fm6gfDvelY55NA6CN4Img4CEdC3tdysm/fmZLZ4STQ3Qmiq98UQh4pGYPWjaQwF5Ho5vGcDscem/990NrSxYmgwZnTH3oSaDwW2Beyffzt9u8NrsvTk666kRZ+BeZ+5sx+ji4ozEUkcSQlnx4EPhPWOmvztIV+xxNBV78pBI4zkenCUZiLiPSXMc70y5Q0GDrK7WoA6NMiD8aYAmPM0jM9RkREIqPXMDfGlAHVgM8Ys22gx4iISOT0GObGmALAZ62tttZWAnWBff06RkREIqu3lnkhsDdkuxrIG8AxIiISQb2FuReoC9n2Bfb19xiMMUXGmK3GmK21tbX9q1JERHrUlwHQ7A7bvoEcY60tt9bmW2vzR42KjtFfEZF40VuY+2jfyvbidKP09xgREYmg3sJ8DU6feFBeYJCTkEHObo8REZHBYazteQEbY0xp4OkRoNJaW2WM8QL7gMnWWl9Xx/TymbXAOwOsOQc4PMD3xir9zIlBP3NiOJOf+WxrbZf91L2GebQxxmy11ua7Xcdg0s+cGPQzJ4ZI/cy6zbeISBxQmIuIxIFYDPNytwtwgX7mxKCfOTFE5GeOuT5zERHpTEvgRrHArKEVANbaEnerEZFoFlPdLIm0zK4xZj5QDxS5XYtEVmCpi3XGmFJjTEKta2SMKUugf9MVHR7ecH5+zIR5Ai6zWw2MABKmRW6M8QZCrd4YszcRVt8M/IxTgGU4V09XuFrQIAr87HH/37iDNcGHtdYXzg+OiW6W0GV2gWpjTJ0xpiCerzQN/oc2g3Qz2ChxI1Bqra00xhQD64C4/wsIdqEZY0qARGioBM0HerzAMM5ssdaujdSHx0rLXMvsJobK4AnaWrva7WIGQ4cGST5Q5lYtgynQtRIMNp+LpQyKQJfKwsBvnusi0bUUK2HupQ/L7EpsC/zmBbT9NhaxVkw0CXQvFQOlQNz+ttlBXuh/7wSxxVpbiHPCrgj3+EishDn0bSleiR/LgdvcLmIwBLrUynHGR7YFBr/jljGmiASbX26t9QV/2wx0tVQS5vGCWAlzH4m9zK7X7QIGU2DhttvCPUAUzQL/2Ctx/pHH+1olhcD6wESGAqAs8JtJIqkmzBkWK2GekMvsBvrVlgOJNCWzFFiVKEHeRYjlAY+4UctgsdYus9YusNYuwDl5lcT7GElXJ6twZ1jMXAHa32V2JfYEpp92nFe/Op4vmAo5Ya/D+Q30kQQ6kRXhjBNsxQn0uP03HRgDWg5swZmZF/ZuppgJcxER6V6sdLOIiEgPFOYiInFAYS4iEgcU5iIicUBhLiISBxTmIiJxQGEuIhIHFOYiInFAYS4iEgf+PyXRMi74pGJBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(T), dAubar, label='assets')\n",
    "plt.plot(np.arange(T), dClbar, label='MPC')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d52d8e84",
   "metadata": {},
   "source": [
    "Save for use in our section 3 plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1b8cb35e",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt('impc_lb_italy.txt', dClbar)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bf81468",
   "metadata": {},
   "source": [
    "# Figure C.1(b): Verify FOSD property required for bound to be valid "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daf93144",
   "metadata": {},
   "source": [
    "We need to verify that for all $t$,  \n",
    "    \n",
    "$$ \\mathbb{P}_{I}\\left(\\frac{z_{i}}{\\mathbb{E}_{I}\\left[z_{i}\\right]}MPS_{0i}>m\\right)>\\mathbb{P}_{I}\\left(\\frac{z_{i}}{\\mathbb{E}_{I}\\left[z_{i}\\right]}MPS_{ti}>m\\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd19b4fa",
   "metadata": {},
   "source": [
    "We will check this in our HA-one calibration. First, we obtain the steady state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "87654d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "with open('solved_params.json', 'r') as f:\n",
    "    params = json.load(f)\n",
    "\n",
    "calib_ha_one, _ = calibration.get_ha_calibrations()\n",
    "ha_one = models_heterogeneous.ha_one\n",
    "ss = ha_one.steady_state({**calib_ha_one, **params['HA-one']})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffdef5c5",
   "metadata": {},
   "source": [
    "Next, we calculate MPCs at each idiosyncratic state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8e323abb",
   "metadata": {},
   "outputs": [],
   "source": [
    "internals = ss.internals['hh']\n",
    "a, c, a_grid, r = internals['a'], internals['c'], internals['a_grid'], ss['r']\n",
    "mpcs = np.empty_like(c)\n",
    "\n",
    "# symmetric differences away from boundaries\n",
    "mpcs[:, 1:-1] = (c[:, 2:] - c[:, 0:-2]) / (a_grid[2:] - a_grid[:-2]) / (1+r)\n",
    "\n",
    "# asymmetric first differences at boundaries\n",
    "mpcs[:, 0]  = (c[:, 1] - c[:, 0]) / (a_grid[1] - a_grid[0]) / (1+r)\n",
    "mpcs[:, -1] = (c[:, -1] - c[:, -2]) / (a_grid[-1] - a_grid[-2]) / (1+r)\n",
    "\n",
    "# special case of constrained\n",
    "mpcs[a == a_grid[0]] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b72475b9",
   "metadata": {},
   "source": [
    "Now, let's evaluate how the distribution moves over time, starting from the income-weighted distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "894f055e",
   "metadata": {},
   "outputs": [],
   "source": [
    "D, e = internals['D'], internals['e_grid']\n",
    "De = D*e[:, np.newaxis]\n",
    "assert np.isclose(De.sum(), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e8d9b3e",
   "metadata": {},
   "source": [
    "Will use some undocumented parts of `het_block.py` to get a `law_of_motion` object that gives the action of `hh` on the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "6cb0ddee",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sequence_jacobian.blocks.support.het_support import CombinedTransition\n",
    "ss_hh = ha_one.extract_ss_dict(ss)\n",
    "law_of_motion = CombinedTransition([ha_one.make_endog_law_of_motion(ss_hh), ha_one.make_exog_law_of_motion(ss_hh)])\n",
    "assert np.allclose(law_of_motion.forward(D), D)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a0e6b9f",
   "metadata": {},
   "source": [
    "Display distribution of MPCs starting with `De`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ae5e0328",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_cdf(x, D, x0=0, *args, **kwargs):\n",
    "    D = D.ravel()\n",
    "    x = x.ravel()\n",
    "    i = np.argsort(x)\n",
    "    x = np.concatenate(([x0], x[i]))\n",
    "    D = np.concatenate(([0], D[i].cumsum()))\n",
    "    plt.plot(x, D, *args, **kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e00fa716",
   "metadata": {},
   "source": [
    "Plot Figure C1(b):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "b202179d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEaCAYAAADZvco2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+YklEQVR4nO2de5wcZZX3fychiRBIOpMLBCGQDiLwsiHMTBAVBLEHd/2IsDhDXkVcRZlxlddlL06Dq6/uu2icyHpZ3YVOXFd33cVkRljlshtmgogL7pqZCeAFEKYJ1wSSzHTIQi5zed4/zlPT1TVd3dUzXX39fT+f/lTVU09VPU91dZ1+znnOOWKMASGEEFIMZpW7AYQQQmoHChVCCCFFg0KFEEJI0aBQIYQQUjQoVAghhBQNChVCCCFFo+xCRUSiItIuIhGf/VnLawERaSzx9XLea5KJiEREJBb29+R3/lJ9X6V+DkltUzahIiLdInIqgBiARI6qURHpCnC+iIh0iYgRkZhPnZjd3y0iUVdZty3vEpGEiAyJyJQ22R95wtYfsPW6C/3R2/N0AxjwlCeyXXcmOOe0bcx3r2d8nTDOXUwKbOdNALoBNIfUlqzPgd0XQYjfV77rh4nrt9olIr05fq9V8UwFxb5reu27JuEW5p593Tn+aJT0u5oWxpiSf6A/1EbXtgEQyVG/C0B7wHMbAF0++3qzXQv64x1ybUdsvZirrB1AwnNcxPbFt+052hnV2z+lHbEAx05pS466MU8/ct7rgOfMdi8Ctb3cn0LbCX3h+j57hXwXQZ8Dz/4Zf18zuX6x++u6pxH39bP1sVqeqQL7HsvR38Zc37fr2M5y9yPX56j8Yqe4WAkcMcYMFnDYeuiDuDHPuWMA+qA337uvFSoAYsaYlGd3I4DJ9hhjUiKShD7wzrEdxpgm90G2XiLL+aaFMaYvYNWm/FUKPmchTLl+SNcpOtNo53Ce/YG/ixphRv11/oE7vxljTNL+1q6C5/ddLc/UdJjmO6MNwAYAHXZZkZRD/eWoFLw4qqUhOwyMOjvsF5AUkfY8524DEIcKCS8tAPp9jlsFIOltj6v+JgDXZTuwkAff6scTdgib8Oxrtf3udpU5KoKEiPS6qjdA71erVWFMHmuvMSIijdnO6fQt2722x/Y66karpnDqRYNc39V2R8XRadvlqOBgywbs9drteq/nGnDVHXKrC+z1Wj1t7HXtz6peCdDOAbve7lHJRHK0M+Ne5Lp+kOfAh1zfl3Mdp++ddrvV72T5rm+/0073szCT/nqIYepvbRD6G3S3IeO7CvrMuJ63Ic9z7PcsTj7z9poDtk7E/oZ6bZn3OyykzzPG9nMI+gc7GvR6M/0dTqufZRj+ZaiVXGUJqDopAlVTDXnqJAB05zl3r10OIFPl0wkVEjHveV31W13b7c61YIekRei3V+XX7j2vu48AWt39hVXp2T4M2HvUbvsVtdsjtq9dzrW89y3fvbbtTLi2HfVE1O/6PtcZcvY5343z/bj2D9jzOe0Y8Ll3A3AN+Z3ru7ZbYVVU9ly+6pVc7URaveDuV69fO32+i3zXz/scFPLbsO3vcm1H3Pdmms9hr6vtQ869n05/s1y/C1NVpwlvWY7vyveZsdvuZy6R71mE/r4HbJ1WV3+GoL+liG2zcR1fUJ89fXKesW6n304/4aOed90L9zPk+x17jpv273C6/ZzRi7LQD9K2iqin3Hg63pilrBM+Lx1Xh52HqAvpl3AE6R9Fa7Zz2GvFkH4hu3+krQBGCuhfNEv5FGHmfEFZfnBuoTIC/fFGPPWm/Ajtsdn6NnnOIPfa59zeY/yu7267t7/e63gFRcx7P7J99/a+dcEjKFx12733wdZtzdLOKT+SLNu9udrpvRe5rh/0OSjkt2HvjVvIdMIlNKb5HLr734XMl1Dg/vq0IeN5dJ1zih00y7Pr+8zY+9vr7UfAZzFbn0a898hVv6A+Z/kO/ARHVpuK7UfCr/05rjWj3+F0+1lym4olm556sswYMygi3nop6M31IwZ9AcAunWF7lzGmw643eK/tDAWNqrGyqbLyXdfNTQA2ZylvtOcJjDGmR0TWwv6zE9U7N5nculivWsGPfPd6pqzF1P46bYva9WFPnVxt7wHQZb+rDmNM3Kp72qG65SFX3SaouqrTbqegatEp361J6/O7RCTunC/LPQ7aznzXb0eBz4El1/e1Edr+RqN2ylXGGD99e9DnMGZVIEGeJ9/+is7cavfU3Yep6ukGZH6HfuR6Zlo8+2DU5hnkWfTiPJ/OeZL2njsEfsb8yPM79tIO6Gw4TxvjUPsKfO71TH+H0+pnSYWK/ZIB7dAgkN0PxXmgPTc+m93DTZv9wBjTZ/WhrcicLhnB1JvcnOe8/bZNrcaYHr9Kkp4Cug8uo79l0ujvaYsvIhI1xsQBxO396Eb6JTotCrjXOY8JgKM6yIafXcsXk2nMdUhAX6YpqNBxSNljAt0nY8wqUXtVO4A+U9gEkmz4Xl9ckz9cRPxOFOT7sr+pjQA6rP0hm73SIe/1rc68yxizym43ILuN0sFpR7b73WE/3vZ72xCF2gpmwhCyTNBBkZ9FSwoI/owVgRZjTIu7QNS+nBCRuDEmZf84e+91M2bW9xRQeD/LYahPIvu8f/eDFocVEJ79uX7wDZ4XYx+Am4wx7hkl2QRTxswvL/acGwBsEs/ccdG55a2uekmfL6APKvHdRs+OLPXctLpGUUnbbnc7G/Icn4tc9zoF/acasds3+Zwj1/W3AJMPvsNVADa6vqMGBB8BAnoPHf0zrIBvANBm74/DZqghc/LH5DXmurH/8HqNMRt8BEqQdrrvRa7rT+c5APL/Npx/qR0m98SRINd3+05E7H7vdx20v1Nwvje3kd+WBxHmub6LHmT3aQvyLHqJIPfznbPP9g9tzgkYQf+s2fdLr7fc9V7LNXlppr/Dgr5bd+NK+oG1W3j0lb1IG4wSyG6XGEAWXTHSviIjcPkT2Ou4jfXdUH3gENL69VbXdXP6wcDqF239Xnu+Vk87fOfvw+qz7acbad2mY4Rtde1vt59eu7/Lc61GV133sSPI1I96z5n3XiNtnDVIG/Am25nn+kNIG8wbXddo93w37batkxMqkNZb+/kYxTBVZ55Aljn7SBuUh9x1fNrp2Gfcn+6g7fTei1zXD/IcZPku8v42XL+PvP4L+a6PqYbqTnsPYtPpr08bGu21HRtmJEudbL+HIN+F+zfa6SnP9iy6r9Nq6zoGc/f35twj74SVbN9xL7JPPIi5zp1A5oQJ9z7ne3H6nDG5wvV+G7EfXxtHnr4HuacFfbfGGIg9sGRYCT1g7PA64DFR6A+9Yn0CrNogajJHRqSCsc9iN4C4sf+U7Wh0AGq/mqkqrGSISJdRdSkpMyIyhPz2z5ql5Oove6O7JL/PiZsu+PiJVBARAE3TtEGQ8tAOVVm6HV8HoarGVLkaNU32lbsBRH0/oGrIVLnbUi7KYVNx9IGRIC9gOwJIVPq/RmNMjzGmrh+mKmSKfc+OVAZNpp2mIrGOaRGrv/edREJKykZTw5EAglBy9RchlYT9Z7kKrimtpnSzemaE868Yqv+m2pVUBBQqhBBCika5nB9nxJIlS8ypp55a7mYQQkhVMTAwsNcYszTMa1SlUDn11FPR3z9dvyVCCKlPROTZsK9RFkM9IYSQ2oRChRBCSNGgUCGEEFI0Qhcq7iQveepNxtEihBBSnYQqVKwj2QhyBz1zgvolAaREs9ERQgipQsIeqSQBLIJGVs2K9ZhPGWOS1hN1OOzUnIQQQsIhVKFiNM5/Kk+1FmQm6MmW84EQQuqanTuBt78d+Pa3y92S3FSCoT6CABkeRaRdRPpFpH/Pnj2laRkhhFQIjz0GPPwwsGlTuVuSm0oQKsDUhDgpbwVjzEZjTLMxpnnp0lAdQgkhhEyTShAqKWSOTCIInmu97unr60NPDwPUEkIqg1IKlYh7w2WM3wy1qzhEazF0dCqVQjxe3BxKHR0diEajiEQiaGqq2PxlhJA6ohR+Kq3Q8NyTfihOxj0Ridg8KYPWl6UTOWaKVTP9/f1IpVJFO19fXx8ikQii0ShisRgaGhrQ11dzspgQUmWEHlDSGNMDTwIhOyNskWs7FEFyww3AI4+EceY0a9YA3/hG/nq9vb0AgMHBQTQ2Ns74ur29vVi1Kp2RORqNIpmk1pAQUl6qMkpxtZFKpTA4OIi2tjZEIpGsdTZs8M8LFY1G0dqaGWwglUqhoSE9vyESiRR1JEQIIdOhpoVKkBFEKYhEIhgeHkZ7u39ggc7OzoLPOzw8nLHtJ7AIIaRUVMLsLzINvCOTVCqFaJQ+o4SQ8lLTI5VKIZlMTto8/F78haq/1q1bh3g8PjnCSSaTiMUY3YYQUl4oVEpAKpVCMpnMafMoVP3V2NiIxsZGxONxLF68GF1dOYNAE0JISaBQKQGNjY0YGCh+8GUKEkJIpUGbCiGEkKJBoUIIIaRoUKgQQggpGhQqhBBCigaFCiGEkKJBoUIIIaRoUKgQQggpGhQqhBBCigadH6uYtra2jO1NmzYxqCQhpKxQqJSIVCqF9evXF90Lft26dZPrFCiEkHJDoVIiip35EQDWrl07JdAkIYSUk5oWKl/6EvDEE+Fe44wzgL/8y/z1ip35MZVKYfv27WhpaQGg+eopYAgh5aamhUqlEEbmR0BHKp2dnejp6UFbWxuGhoaYU4UQUlbEGFPuNhRMc3Oz6e/vL3czCqKpqSmUSMUOLS0taGtry5ldkhBSvfzkJ8DllwOrVwOPPjq9c4jIgDGmubgty4RTimuEaDTKUQohpOxQ/VUCwsj8uGHDhimJvZj5kRBSbihUSkBYmR/b2tqwdu1aRCIRJBKJGbaSEEJmDoVKCQgj82MsFuPIhBBScdCmQgghpGhQqBBCCCkaFCqEEEKKBoUKIYSQokGhQgghpGhQqBBCCCkaJREqIhITkZzRDm0dzpElhJAqJnShIiIJAEkAKRHJ6qwhIt22zrBfHTKVVCqFeDyOeDw+ZV9fXx96enrK0CpCSD0TqlCxI4+UMSZpjOmDCo1YljrDts6grRMJs13lwBEAxWJwcBCLFi3Cxo0bp+zr6OhANBpFJBJBU1NT0a5JCCH5CHuk0gJgyLWdBOANfpUE0C4irSISBdBrjEmF3K6SU+wkXdFoFCMjI1MySfb19SESiSAajSIWi6GhoQF9fX1Fuy4hhOQibKESATDs2k7ZskmMMUkAbQC6AQwAGMx2IhFpF5F+Eenfs2dPGG0NFXeSrmIQiUSy5mbp7e3FqlWrJredQJaEEFIKShH7q8GznXJv2NHJOmOMiEg7gF4RabKqsEmMMRsBbAQ0n0qQC2/aBIT9Po1Ggeuuy10nrCRdftdqaEjf8kgkUvQ0xoQQ4kfYQiWFzJFJBKruchMHkABUcIhIE4Bm+IxYqpFIJILh4eGcCbQKjVKci+Hh4YxtP0FGCCHFJmyhshlAFwDnb3jUGuwhIjG7nkKmnaUBQFGMAPlGELWId2SSSqWYvIsQUjJCFSrGmEERGRSRLgD7oKMS2Nld3SKy0hgTF5GEiHRCBUzC2llqhjCSdLlxC5F169YhHo9PjnySySRD5BNCSkboNhVjzJR5tHZ21yLXdkfY7SgnYSTpAoCenh4kEgmkUin09PSgtbUVjY2NaGxsRDwex+LFi6fMDiOEkDARYwLZvKceKHKJMeb+IrcnEM3Nzaa/v78clyaEkLLwk58Al18OrF4NPPro9M4hIgPGmObitiyTQCMVEVkJ4DakbR8CYCWA2SG1ixBCSBUSVP0VB/AMgBvhmRJMCCGEOBRiU+k0xrwaWksIIYRUPUE96rsAZDhZiMhfFL85hBBCqplCwrRsEJFx+5mAChpCCCFkkkBCxRjzDIC4MWa2/cwC8Ilwm0YIIaTaCDxSMcZ81bO9qfjNIYQQUs0EEioiskBE7hORCav+2hx2wwghhFQfQWd/fQdALzREPQA0i8jHjTHfCadZJAipVArr168HgAzP+ba2tox6mzZtYlBJQkhJCKr+2m6M+aoxZr/9bENmnhSSh1JmfgQ0BpjzoUAhhJSKoCOVkSxl6wDcUcS21DRhZX7csmULhoaGMvatXbs2cP4VQggpJoGFirWjJKE5UWKwOVAqmR//GHjppXCvceKJGo8nH+7Mj42NjTO+rt/oI5VKYfv27WhpaQGg+eopYKqLkRFg717AGP0A2ddz7QtaL8i+iQlgdBQ4ckQ/7nW/snzbuepMTEy9J34hCmux3K/uwYPZyyuNQELFGPMjEUlCRyf7odOJt4fZsFqilJkfAR2pdHZ2oqenB21tbRgaGmJOlQrm9deBxx4DduwANmwAdu4sd4tmxqxZwNy56c+cObm358/P3DfbJ6KgSH2X79oF3HUXsGxZ9vqVQk6hIiJrjDGP2E0D4Ieu3V0A/jikdhWFICOIUlDKzI+RSGTyXK2trYjFYujr68t5bVIYO3YA990HjI0Br7wCPPCA/st2/8Mv5LNrV/rf+YIFwMUXA9deqy9nkfTLJd96WPtEgguIXEKBzIyHHgLuuQc4+uhytyQ3vkLFRib+DnSm10Joet8BaIRiQKMUV7RQITrK4SiluNx4owoVh/POA04/feqLOOjnpJOAc88FGhuBk0/2/+dKSDXgK1SsF32zXd8vIi121hcAQETeVYL21QSlzPy4YcOGKaMeZn4sLgcPAhdeCGzbpv/KZxUS7IiQGidoPpU10HTAbhYWvTU1SqkzP7a1tWHt2rWIRCJIJCp+PkXVMToKHHusqnoIIZnks6ksgCbmarbbKbsrAuAmcEpxIBobGzEwMFD087a2tk4ZwcRiMY5MQmZ0lAKFED/yjVRWAei26w3IdHjsCaVFhFQ4Y2MUKoT4kVOoGGN2ADgNUBuK26ZCSL3CkQoh/gQNKLkZOtvLXcYkXaQuoVAhxJ+g81YkS/DIVcVuDCHVAIUKIf4EFSrJLGXNxWwIIdXC6ChwVNAAR4TUGUF/GputCqwXOvOrAzTUkzqFI5XCGR9X/56DBzUsTa718fHMOGTO0i9W2Uz2Z6vjV7/U1/LuDzuOYbEIGvtrh4jcCBUmUQA3GmN+FGrLCKlQqk2oGKMv6sOH9XPkSHrprB8+7B/IcGIi/cLPJxCyrb/+ul6nHHhDz3iXfmFqZrI/Wx2/+n7tyXaOAwd0hFzVsb/cWA/7G51tT1wwQuqGMIXK0BDw5S8DjzyiL3Pvv1tnGeSfsvuYMJg7V+NQHXOMLp31hobMcu9+Z5mt/Oij06rFmQqBWuPxx4F3vANYsqTcLclNUI/6hQC+gkw7ShTA4jAaRUglE4afyvg4cPfdgBNY4QMfAObNC/YSzbXubDuBH+fN04+z7l76hZuZNQt4wxumCgkGjiTZCDpSuQkaSHIL1JbSCHWGJBVKX18fUqkUc6mEgHekMjqatgNk+0xMqCBycoaMjelydFRVRP/+70BvL7BnD3D88cCHPwx8/OPl6x8hM6GQOSyd0GnE+2x+lVtDalNN4uSTd+eSn+n52tra0N/fj4aGBiQSicnwLB0dHYjH40gmk2hqagolREy94tgndu0Ctm4FHn4Y2LIle2KpoMyZo2qNyy8HWloYoJJUN4FnfwFoN8bcIiJbROSX0OyPgRCRGICIMcZ3xpiIRAC0A0jmqlcIDz8M7POGwSwyixcDb3tb/nrFTie8ZcsWxONxxGIxbNiwAS0tLTDGoK+vD5FIZDLkfUNDA/r6+hgPrEiMjuryrrs0vwUAnHMOEIulhYHX4DprVlr9dNRRuj5nTnr9jDMqX09OSFAKmf2Vsptx6Cywq4IcKyIJaEKvqIgMGGOastSJANgE4DpjTMq7vxYodjrhWCw2GUa/s7MT8Xh88jqrVqX9Up2Q+6Q4OELl+OOBO+5Q+8Ib31i7xmFCCiWoof5KWAdI7yywPMfFAKSMMUkASREZFpGYMabPU7UbQFuxBUqQEUQpCCOdsDsvS19f3+T+VCqFhoa0uSsSiRR1hFTvOELl2GM1MRchJJOg6q/PAsh4q4nIqcaYnXmOawEw5NpOQmeNuc/j/G2/SkRaAGzOpv4SkXaoegwrVqwI2OzKIOx0wolEAps2bZrcHh4eztjvJ8hI4ThCpZr8VAgpJUFNgnEAcRG5RETW2KRd8QDHRZAZLj9ly9zEoIKmz55zkx3hZGCM2WiMaTbGNC9dujRgs2ufeDyOTZs2TQoO78gklUoxnXARcYTK3LnlbQchlUrQkUo/gASAta6yoDnqvVOPU57tVQB6rIoMItIHHeF4VWRVS1jphOPxOG666aaMkci6desQj8cnRz7JZJJG+iIyNqZLChVCshPUUL8fNq+KQ8Ac9SlkjkwimBqcMuXZTmJq6uKqJox0wh0dHdi4cWOGMOrs7ERXVxcaGxsRj8exePHiok1hJoozUpk3r7ztIKRS8RUqIrIVQIcxZqdVdzUYY+53Vdke4PyboTO/nDdf1DHSuwz2Th2HRujsspohjHTCiUTCN/88BUl4UP2VHWPSMcQOHZr68St39h08mK7jDk/jF6bGXVZI3XxlMz2HNyxOkLKgbTh4EFVBrpHKBgCtIjIIfdEvFhHHPtIAoA151F/GmEERGRSRLujoIw5MTiHuFpGVtk6viHRCRy0JRxVGSKVRipHKoUPA88/7x/dyvPSzlfvFBxsbmxpE0r0+OurvwDkxkVsw5AtImY3ZszX0y7x5unQ+8+dPTSuQLexMtmWhZTM9PlvbgpRN93r79wPPPQe86U2oaHyFijFmmxUiN0EN6REAjo9JAwLaVIwxUwz6durwIte2v0GBkArCibYbllAZHAS+9jV9gZSCWbN01JUv9pf75X/00cCiRbp0l8+bp2XZhIW7zhveoIKDvj2FsXs3cO+96htVyQTJUX8VAIjIuXYbdjuITYWQmuKZZ3S5aFHueoVy4ABw333A976nL9143D+gZK5ovX77jjoqM3ik82GyMVJsCgl9v8Ozva34zSGksvnlL3W5cmW6bPduVR85ailvMMmJibSKyJ3TxIkj9vDDwK9/rec66yyNAVYpjruEFAr/pxBSAAMD+g//2WeBf/s3zXsy0zkYy5YBV10FXHBBprAipBqhUCEkIMYAO3aoOqmvD/jtb7X8ggt0ZOEEjwR06aieZs9O5zFxPnPmpOsuWMDIxKR2oFAhJCDPPQeMjKihdMUK4Ac/UMM0BQIhaQL/HETkVBG5JMzGEFLJfOc7upw9WwXL/PkUKIR4CRql+DPQqcVDsKFaROQvjDG3hNg2QkrKU08B3/oW8PrrU/elUsCPfgS8973Af/4nnR8J8SOo+muxMaZBRK5zl4XRIELCwBjg1VfTHtsvvKAfJxXwk08CTiCCxT5PdkcH8O53Az//OZNqzQRj0s6Yo6NpJ0wn3bLbCTObd7l7Pd/+sI4Ncp6g1wh6nlL5Ls2UoELFTqSEu6sx6OiFkIrl5puBv/s7YO/edDBIP/7wD4Fvfxs48UT/Ov/3/6pvx/LlxW1nKXGmMo+O6mdsLD0l2q++88L3CoBCtt3lpHCc+7ZwYXnbkY+gQkVE5DYAK22IlXWooSjCpLp57TU1oE9M6Aty505NI/3LXwK33AKsXg189KM6AnG8kU88ETj5ZFVjOTaSU07Jf63+fnV8DMuj/rXXdGbZzp2ZoVi8S791v6VXiBSbOXPSs9rczpXz5+v9cra9+7Ntu9Myu5dB1/32Tfc8hV4njPOIqFr2X/+18qedB41S/CObTrgNGq34K8aYH4XZMELysXev/vNdswbYsyd7nXPOUTvJhRfO/HojI8DQEBCJhJOk6+WXta2HDgGnnabX8Itz5ffi9fOynz1bz3fUUbp0PkcdlS7zC5si4i8A5s1jyJVSMWeOfo+VTlBD/anWg55e9KSsjI0BV1wBbNumL1+Hj30MeOtb9WW7YIGqp04/vbi2j7vu0uVxxxUW3sTxrB8b0xHD2Fh6/cgR4De/AR5/XINIAsDb366qOEKqkaA/jTYRMZztRcrJE0/oP/l77gE++EHgzDNVaMydC6xbp6qWmeCETHn99akhV4aHga98RV/4u3er0fSJJ4BkEviv/9Jj3eFZvJ9ciKgq7j3vAc47Dzj22Jn1g5ByElSoJAEM2tlfBsAWY8yr4TWLkDTPPw/8xV8Ad9+tL/zf/311PCymymX3br3G9hxZgs44QyMIX3yxjpRefFHLjz8eePObswd8dAd0nD07c+msn3QS0ODNj0pIlRLYpmJXNwGAiNwqIg3GmHWhtYzULb/7HfAHf6CjgddfTycnuuwy4OtfV0NlIQJlfFwDNh44oKOG0VE14j/yiJ774EHglVe07s03qx0mm2A46SQ91/i4jpD+z/9Rm8KyZXSCJMQhqE3lEmie+nZoVsb9ANaH2C5SR+zYAfz5n+uMrf37NVgjALS06At+/nygsRF43/sKP/cddwDr16uPips5c/TcJ56ooVYWLAA+8hHg1FNzn+/QIRVM8+cHmy1WLzi+J26bUa6PX73x8ex+Kvn8OopZN8g5CqlbrLY525U+KSKo+utS6BTiBICrvGHwCSmUAwd0eu7XvqZqrYYGnaG1YIHOw1+5EvjTP53ZD+i554CbblLB8dd/rWoqJ+jjqlXTs8E4/hxhetQfOKCz2Yr5QpyYSI+y/D65/FRyCQvn+OngVgU66kDvzDa/db9pue5RY766hZy3lHX91mfP1rhzlUxQoTIEYJExZtKnU0QW0K5CCsEY4IEHVJh88YvpcChtbcAnP6m2ipkwMaHhVO6/H3joIc2SBwA33KB2mGLgeHyHJVT++7+Bxx4rLDXvdJg9e+onlwB3Xvpz56pPj1sQeO1EQffluyapTnIKFRFZY4x5xBizSUTWSOYT0IEA6YQJGR8HvvAFoLtb7SWAjko+8QmguVmn/k6XiQmdfvvUU1P/LS9cqJkUzzpr+uf3UuyRijE6SeCFFzSrZCqls78uvXTqS3e6/9yd0ZkjPGj/IWHiK1REZCWA7wBoFpGFAAYBDABwHt2VoFAhOXjmGeDOO9Ve4nDDDcCnP622i+n+S/3Nb1RlNjCgNo4nnwQuukg95+fOBZYu1e2FC4vvLHbkiAoyERUAe/bo9Z3px0D2paN+cjJBOutu1dEJJ6iwXbFi5tOjCSkXvkLFGPMMgGa7vl9EWtwphJmjnuTilVfUK3xiQmdKNTbqNOClS2d23hdfVJ+U0VF98b7lLWp/+cxndHZW2OzapWq73/0O2LJFy+bO1T7m8m53RgjOiMFZnzVLjz3llHC89AkpNUFnf60BsM9TXOFhzUg52bdPBcrNNwM33jjzEYMxGuzx9tv1Jf3jH6swCSsGlx9OLvk3vhG45BIVBMuXMxQ+IQ75bCoLAERhRyw2/hcARKARiu8IsW2kinFCqJx99swEyvg4MDioPiEjI8D556s3/RlnFKedhfLrX6vN47jjdCRGCMkk30hlFYBuu94AYNi1ryeUFpGawBEqb3hD4ce+/DLw2c+qvWLnTs27AQB/+ZfANdeUb8bQgQPanvnzqaoixI+cQsX6o5wGqA3FbVMhJBeFCBVj1AHyuefU+P2LX2h2xWgUuPpq9Vg//XSNu1VO7r9f21poQElC6onAPw2vQHGmGxe9RaQmcELRBxEq//APwFe/mt5esAB429uA224rrc3kyBEVaI7XvDOja2JCPfI3b9Yw+zt2cKRCiB9BDfULAXwF1rYCnVa8EkwpTLIwOqrOjEDu0POvvQa0t6sz5LJlwHe/q/UjkdKpuJyshkNDGgH5hRf86559tnr5f/jDFCqE+BF0pHITVJBsgdpSGqE2FkIyGBvT8PH79mlO91WrptaZmAB+9CPgpz9VgXLiiTqr64QTwmnTyIgGj3ztNb32kSNqcH/8cQ0m6fiTLFmiTprRaNph0B1Qcv58HbEYQ/UXIX4U8tPohBru99lMkLeG1CZSpRw+rEEgf/5z3Y7FMvcbo4LkttuARx/VabjveAewaVN4bdqyRf1jvGFPli/XsDALFqiK7thj1fHQSTfshxMjiyMVQrITVKhsBtBujLlFRLaIyC8BtAS9iIjEAESMMTlnjIlIAkBvvnqk8vjNb9TeMDamI47vfEcFjENvL7Bxo8a1AoAPfQj43OfCVXPt2AH88z8Db3oTcP31mqPecT48+ujpXXtsTEc7FCqEZCdoPpUdLh+VODTu17VBjrWCogtAVEQGjDFNPvViAGIAeoOcl1QO99wDvPe9uv5nf6YZEp2XrjHAffdpaBZAsxt+/vPFTUp15IhO9z1yBHjpJVV1PfqohomZP19D2kejxbnW2BjVX4TkIlfsr1OhTo7usjV2dTM0t8qDuU5uBUXKGJMEkBSRYRGJGWP6slRvhMYXI1XGv/6rLrdu1dGJMwJ4/nkNn7LDJkrYskVzmBQLxz7S3q52EzfHHaeBJD/3OV0vFk5eFsbmIiQ7uf5vtUFHJCmf/UECSrZAw+Y7JKEe+hmISCt0AsDaHNcjFci3vqVC5eqrNbKuw969aZvKhRcCn/pU8QTKtm06+nn88bStpKUF+L3fUzvN8uXFG5l4GRpSoXniieGcn5BqJ5dQ6QOw0Z1DxU3AgJIRZHrhp+AZ/ViixpgeyaHkFpF26OgIKyo9S00d8cADuvz61zPLd+3S5Wc+A3z0o8WJFjw6qvlGvvENHSmcdRbQ1KS+LO96V2lGD08/DSxaNL1IAYTUA7miFOfL7hg0oKRXe55yb1hhsTHfSYwxG516zc3NIacwIkF59VWNx+WNPvygVYw2Nc1coExMqArtu99Vr/ulS4Fbby19MMmxMb3+kiXFD6lPSK0Q1Pnxaai6y0u+n1YKmSOTCFQF5qYFQIcdpUQBNIpIwhizIUjbSPkwRo3ijY2Z5YcPA9//vjoxvulN0zvvCy+or8vzz+sU5ccfV9XW9dcD551XeoECaLj7sTEVKkx0RUh2gs5h6QGw3lGFWRvIogDHbYbO/HIERNQx0jsGe2NMm1NZRLoBbOaU4urgjjvUduKeyXX77ZoqGAC+/GX1/ygEYzT51l/9Vbps1izg5JOBDRsKP1+h137xRTX+G5P5ef11oKdHBeXSpRypEOJH0CnFN3q2e0RkK4CcbmvGmEERGRSRLmg+ljgAiEgEQLeIrDTGpGxZO3RKcUREksYYzgSrYO64A2ht1fVvfENfvLfeqoZ7QKMMB80Lv2sXkEjocv9+9XwHNGPkmjWawTHssC179qiX/9NP+9dpaNAJCVu3cqRCiB9B1V+neopWQWdq5cUYE89SloJnpOO2mZDK5sgRoM2OL++8U+N2/c3fqGf86tXA+vXBc408/bQGk3zpJVVrrVmjo5LTT59Z7no3Bw7orK1Dh9Q+MzamIeyff15VdUeOaNkb3gBccYU6b7rDszifE04A/ud/9JwcqRCSnaDqrz4ABun89COwow5Sf2zbpi/nf/xHfQnffrsKlFNO0YjDCxbkP8ehQ5p8a/163f70pzM98IvFfffpx8vChTrt+Oij1VZzzDFAc3P+tu/dq0t61BOSnaBCpcsYE2KEJlItGKOCBNBRxX/9l9pQVq7UUcvRR+c+/sgR4KGHgL//+3TOlVtuAd785uK2c2JCw9jfdx+wYgVw5ZUqMJxAkcccMz2V2ssv63LZsuK2l5BaIahNZYpAEZEFxphXi98kUsns3q2C4eqrVah86lNavn59boGyaxfws5+p3eLQIX3Rt7ZqFOOZuh29+KKOIMbGtH0vvKDqrdFRjfd1xRXASSfN7BoOu3bpOZmTnpDsBI5gZPPVu+POrrMfUkdcd50u/+iP1FbR1wdccAFw7rn+x+zaBdx4IzA8rGqnj35UszguDOrplIXxceCppzR0/b/8y9T9p50GHH88cNllxYvTNTGhI5UzzijO+QipRYIa6q+DTg3uB8Oo1C333KOfSAR45zvVWx7QLI3ZeO014J/+Cbj3Xt3+wAeA979/Zj4m4+PA976nfituLr1U/WXmzFHVVhg2j927dTS0fHnxz01IrRD0P1wLgCZjzDNhNoZUNoN2kvdPf6rRiLduVUHxsY9Nrfvqqxre3hg14F9++cwN8U89pbnrH39cPfVPPlmN7XPmqENiMaYd79unajMg7aMCqMrvwQdVYDHuFyH+BBUqCahH/aRQEZFLjDH3h9IqUnFMTAA336xTaXfvVoFy/vkaxj4br7yiL+SrrgLWrZu+DcIYPdedd+q04KOOAt76Vh3xFIuJCSCV0oyVL73kX2/+fA3xXw5vfkKqhaBCpR/AgIiMQANELgZz1NcVP/uZ/lv/0z8FvvY1Lft//8/fX+ORR3TZ1FS4QNm/H+ju1tHO3r16XUDVbO99b2Hnm5hQobR7t45AjNGyPXvSoxJnNDJvnl4jEkn7pgDpZUMDBQoh+Qg8pRia6ySBtFBpDatRpPL4xS90efAg8Oyzak855ZTsdV9/Hdi8Wf/Zn3xy/nMbo7nq9+1TO8yOHTpDbNky4C1v0Xwoy5cDZ55ZeLsffFBjdgGZjoyLFmlcsrlzVTDOnavGfUYfJmRmFDIv5uPuKcQisj2E9pAKxBjg7/5OX7z3368v3muuyV73xz/WVMKA2l3yJcg6ckTzsfz617p99NEqTE46SacCTyccyuHDmrTr0UdVAJ5yiuaj5yiDkPApZKTSZWN4OXQAuKn4TSKVxr33qq3hzDN19tXtt099QRujDo3/8R+6/clPqu3DD2P0hf/tb+v26tUqqKZjbH/oIRUijm3EcaoEdDRSrqjGhNQjQYXKMNQnJYZ0qJZFoFCpCxzflIMH1aZx1lmZ+41RO8sDDwBr12p4er8c9KOjOnvrzjvVzwXQmWEXXBBcoOzerc6Ne/aoAHFSCZ94okYQnjVLRydLl6qjIiGkdAT1qN9vIwpPZoEUkSLOvyGVyh13qPPiypU6dTfb9OEf/EAFypo16uSYzZD+0ksa0v7hh9Oe7pddpt70K7Nl6smCMSpQ7r5b1+fNU8GxYIGOioLEHCOEhEtgm0qWtMLMvlgH/PEf6/K449TovmZN5v4HHgC2bNHRyxe+MNV7fWJCbRtbtqSFyUUXaT75fPYW9zl27tQZZXv3qk3nD/9QJwIwBD0hlUXYmR9JFbNtm07HXbpU1UwdHVPrOGmDb7ghU6AYo86S//ZvqjYDNLTL//pf+QXB+Djw5JN63IEDOsr5n//RkdJb36qqraACiRBSWsLO/EiqlPHxtKrrmGP040269cILwPbtagj3hi4ZHlaDPqBxvt7yltye6Pv2qQA5ckQFyq5dWj5/fnpG2MUXFy+OFyEkHELN/Eiql2uv1dlZra3AY48BZ5+dGU9r1y51fgSA97xn6vFOMqtrr51q2HdjDPDf/63XcLNggaq4OGuLkOoi9MyPpPq49VYNBLlqlaq/AOC229L7jVH7ya5d6ovS1DT1HDt26HL+/OzXGB/XiL93363bS5boDLB58zQPPTMrElKdMPMjyWDXLvUxAVR99aEPqerJLRzica338Y9rtGI3hw8DTzyhgR+XL5+q8nr1VSCZVKHjBG485xydikyjOyHVDzM/kkkmJtKqqkQiHePrc59L1xkYUD+TY46ZqvZ67TXghz9Mh6X/gz9Iq8ySSS1/8UXdXrJEr7V4sU4EIITUBr5CRURs9nAkvAJFRFYCONcYc0eYjSOl5d3vVo/0s89W7/m/+Rud/vvud+v+V17R1MEA8PWvZ9pYXnpJQ7kcPqxJrC6/XAXH8LDOEHPUaKeeqjPAjj+eRndCapFcP+smY8yl2XYYY54RkU4AFCo1wuc+p1kclyzRsPZXXqnlTiIuQH1NAKC9PVOttXNnOtzKRz6iQmlsTKcT79mj5SKae+XYY0PuCCGkrOQSKsmStYKUlTvvBL70JVVFvfgi8Cd/oqFPNmzQ2FlAOofKypUaqsVhzx4docyZA3ziE+nIxa++qvtOPlmnEy9aVJwkWoSQyiaXUBnJc2y0mA0h5WFsTEclCxYAv/qVGtjvv19nfr3vfel6W7fq8tprM4XD1q06G2zdusxQ+IcP63L1av84YISQ2iPXfBtf50ZrU+H/zhrgm9/U5Wc/q8mpnLAsiURaeBw8CPzkJ5q61x2m5cgRDZ3iLQfU9wRgfhJC6o1cQqVbRG71FlqflfsA3DblCFJVjIykZ3Z9+tOa1REArr46M7nWpk0qQNzBJJ0Q+IAmt3LYs0dD5b/yipYzSjAh9YWv+ssYs01EVonIPqifyjBU5RUD0MOZX9XPv/+7xvS6+26d8vvTn6qqyj2F+He/A3p71TC/enW6fNMm4OmngQsvBC610zkefVTjfY2OqsH/3HNL2x9CSPnJOanTGLNRRJLQ1MGroMb7S40x20rROBIuX/2qBmaMxXR0Mns2cNddmU6IX/qSLq+/Pl12++0qUM47T+0uhw/rccPDKpTe/vapscAIIfVBXk8BY0wfdKRCaoiXX9aRxcc/rvaTX/1K1V9LlqTrPPSQCorlyzVUPaDxwAYGdP2yy9RI/9Ofar3TTwfe8Q56xhNSz9D9rE75279VgXD11ep3cu65unSYmND0wIAa8QEVHN/6lq7H4xo9+KGHgOeeU9XY+eeXtg+EkMqjJP8pRSRmw+WTCsBJ/3vOOcA//qOWXX995gijq0t9TT78YfWCB9L556+4Ih1a5eWXNQgkBQohBCiBUBGRBNQWkxKRgSz7IyLSKyIjIjIkIrGw21TvfP/7aqB/5zuBX/xClxdckN5/112a9nflSg0/D2gYlsFB4I1vVJsJoHaVvXsz/VMIIfVNqELFCoiUMSZpbTPDWYTGVdCAlYsAJAD0htkmAtx8sy7Hx3W5YUN63yuvABs36vr69en4XD//uS6vvFL9Vx58UJ0klyxRj3lCCAHCH6m0ABhybScx1RO/zwocGGM2gITKU08BQ0MaamXrVk3xu2CB7jMG6OzU9S98IR3ufvduzfB45pk6Kunv1/D2J56o5zn66PL0hRBSeYQtVCJQ/xaHlC2bxBgzGWPMjmJ6sp1IRNpFpF9E+vc4UQpJwfzwh7pcvlxtKB/9aHrfN7+paX3XrAGam7VsfBy45RZdv+IKNdYPDur2JZcAc+eWquWEkGqgFIZ6b+SnVI66HQCuy7bDGLPRGNNsjGleygQc0+ab31Sh8bOfaapgx5/kySeBbdtUSHz+8+n6TmTi1avVIH+HdXm94grNqUIIIW7CFiopZI5MIvCJfiwiXQCuM8akQm5T3fL00zoSef113f6jP9LloUMqbObOBb773fToY9++tE/KNddotsaJCc3SuGxZ6dtPCKl8whYqm6F2FYeoYz9xG+ytQFlPgRIu//RPujxyREcrTsyur3wFeP554IYbgIULtcyYtMH+6qt1ltfjj6vHvDd4JCGEOIQqVIwxgwAGRaTLJvWKAzqNGBqwMmKnHHcCGBERYz9dYbarXrn7bp0SPHcu8Od/rmX9/Toa+f3f1zheDt/7no5ULr5YHSN37NDyiy9mXhRCiD+he9QbY+JZylJIh9bvsB8SIr/6lQqGE07QkCvnnach7b/6Vd3/oQ+l6951F/Cb3wArVmie+eee0+yOq1ZlhnEhhBAvjNJUJzi+KPPmqZc8ANxzj9pXvvjFtNrrkUfUiH/MMcCnPqX7HU96x+mREEL8oFCpA8bHdWbXokWaHviyy4D9+3VmV1OTfgD1X/nBD3T9hhs0avEvf6nb55/PhFuEkPxQqNQBDzwA7NqleeRbW9Um8v3va8j66+wE7ldfTRvmP/YxNcjv2KGC5rTTMnOpEEKIHxQqdYCTE2X+fI3l9eSTmnirpUUN96OjGlhyfFynDp95pibn2r5dw7S89a3lbT8hpHpg6PsaZ3xc850AOgJZtiwdyv6aa3Tq8Le/Dbz4IvD+92vk4t/+FvjP/9Q6V17JMCyEkOBQqNQ4zlTgSERHHI89puHq3/MeNc5v3aoC5cILdf/OnWmB8sEPAsceW66WE0KqEQqVGseJ27ViBXDRRTqj65hjdAbYnj1qwF+wQANDPvVUelTT2kqBQggpHNpUahhjgJ/8RNc//GF1dHzxxXRk4R/+UI33N9ygoxdHoLz3vWqoJ4SQQqFQqWF++1t1cGxoUHvKP/+zZmz8wAd0RtizzwLvex+QSqkvypw5wFVXaUh7QgiZDhQqNcz11+vy4otVwDz7rKq19u0D7r0XOOMMnS78H/+hPijvf7/aXgghZLpQqNQoxqhnPKCzvf7hHzRES0tLOrBkS4vGAwOASy9NJ+sihJDpQqFSo/T2qmBZskQdH195RX1Utm9X+8k736n+KmNj6mHPmF6EkGJAoVKj/P3f6/Lqq9PZHteuVcN9Q4PaTXbu1BAtTqIuQgiZKRQqNciRI8CPf6xTh1ev1jhfra3AL36hI5OLL9YRy0knqbMjIYQUCwqVGmTTJl2uWKHOjoA6Nz74oI5Mfvc7zU9/0UXMjUIIKS4UKjWII1Te8x4NCPnBD2pE4nnzdPRy5IjaVObPL287CSG1B4VKjfHUU8Cjj6rd5LXXtGzRIk0H3NysRvrTT9eEW4QQUmwoVGqMm2/W5Vveot7zb3ubqr3GxjS8/cKFwAUXlLeNhJDahUKlhti9G/iXf1HB4XjFn3CCqrve9CYNcf+ud2k4e0IICQMKlRpi/XoNdX/BBarm+r3fA554Qm0pc+YAjY30RyGEhAuFSo1w6BDwt3+rxncnGOT8+Vq+dKlun3VW+dpHCKkPKFRqhG99S5fnnAOMjKi6a2REBcsxx2jkYeaYJ4SEDYVKjfDXf63LaFSXIipU3vhGFTSMPEwIKQUUKjXAnXcCBw4Ap5yiYexXrwYOH9bwK4sWaXgWQggpBRQqNUBnpy5PO03tKQcPArNnq0f9pZeq9zwhhJQCvm6qnIEB4OmndVbXvHkqVF59VW0q55/P2V6EkNJCoVLlfOQjujztNE2wlUoBb34zcPbZqgYjhJBSQqFSxfT1Ab/+NXDccfoBgJNPVn+Uiy4qb9sIIfUJhUqVMjEBXHONrq9cCRx7rPqjvOMdGiyS0YcJIeWgJEJFRGIi0jrTOiTN5z+vYVkWLFC117JlmsHxne+kYZ4QUj5Cf/2ISAJAEkBKRAamW4ek+e53gS9/WUcjJ52k9pRrrgEuuURnfRFCSLkINbSgiMQApIwxSQBJERkWkZgxpq+QOkQ5eFAzON57r26fdJLaTj75STXME0JIuQk7Xm0LgCHXdhJAdBp1itOYFg0DnwtjMpfToZBjC6k7MaHLWbPUS/6664D//b/VwZEQQiqBsIVKBMCwaztlywqtAxFpB9AOACtWrJhWY0ZGNK9IsQlqFM9Wr5Bj584FzjhDc6ZcdJEa5wkhpJIoRWaNBs92ajp1jDEbAWwEgObm5mmNI/r7p3MUIYSQoIRtqE8hc9QRgaq3Cq1DCCGkCghbqGyG2kwcoo4B3hroc9YhhBBSXYSq/jLGDIrIoIh0AdgHIA4AIhIB0C0iK/3qEEIIqT5Ct6kYY6YICWNMCsCiXHUIIYRUH/S9JoQQUjQoVAghhBQNChVCCCFFg0KFEEJI0RAzk3gkZUJE9gB4dpqHLwGwt4jNqQbY5/qAfa4PZtLnU4wxS4vZGC9VKVRmgoj0G2Oay92OUsI+1wfsc31Q6X2m+osQQkjRoFAhhBBSNOpRqGwsdwPKAPtcH7DP9UFF97nubCqEEELCox5HKoQQQkKCQoUQQioIEYmISJcNsuvdFxOR1qDl5aCuhEol3fhyUOv9F5GEt3+13Gf78umssz7HXGkzvOVV32cRaQQwApvl1rMvAc01lRKRgXzl5aJuhEql3fgwsC+ZXhEZEZEh94+v1vtv+xrzlNVsn236iE0ANhpjelzltdznbmjfhiv5pTpDktAI7hmR2+3znTLGJG2+qWGXgJ1SXvpmp6kLoVKJNz4krgLQZYxZBCABoBeom/43Ahh0Nuqgz90ArrNpJADUdp9tP4Zt3wahfYvUWp+NMSn3d+qiBcCQazsJIJqjvGzUhVBBBd74kOhzsmYaYza4ymu6/1bt4fxbT9llzfbZqkgA4CoR6XapfWq2z9C+tItIq4hEAfTal28t99lNBMCwaztly/zKy0boSboqhAgq7MaHgTEm6azbf2vOizaC2u5/1BjTIyLusghqt88x6Iuzz34GRCSFGu6zMSYpIm3QEVoKQJvdFUGN9jkLDZ7tVJ7yslAvQgWosBtfAjoAXOfarsn+i0g7/J3BarLPAFYB6HH+RIhIH/QfO1Cjfbajk3XGGLHfea+INNndNdlnDylkCssIdFS2yqe8bNSL+iuFCrvxYWKnIrr17SnUbv9bAGyzBtoYgISIdKK2+5zybCcB7ENt9zkOtRPCGLMR+keiGbXd54hrfTPSfxwAHZ335SgvG/UiVCruxoeFFSjrPca+mu2/MabNGNNkjGmCqoLi1p5Us32G9q3Rtd0IVXXWcp9TyLSVNEC/75rrs7WRdQCYnCZtJycMWv+VTtjZYX7l5aRuwrS4HIn2QQ3ag7nqVyN2aqV3fvsGY0y81vtvVSJdAPqhgmWwlvtsXyCAvmyHnWnFNd7nBNQonwLgzPaq6T5XI3UjVAghhIRPvai/CCGElAAKFUIIIUWDQoUQQkjRoFAhhBBSNChUCCGEFA0KFVJX2NhRQzbibbb9XTbKc7uNAjskIgM2xHzCrje66sdseZdd9mbLg+G5RqNdtttPq2tfq41AnFGXkGqhnsK0EAIbI6wDQKuIRLJEhG2H+n1sBCZDoKScAJ1WYGwC0GR9Y5qMMR3uE+QSKtbXIm4FScQYs8EKqj7blqinTUkRSXivQUilwpEKqUeSUAe6m9yFVkj0e+qmPNvbbd0YNM1Atpd9IttF7fkHrNBYi3Q4kWEADZ5oywA0FDqAoVpIQEXqAwoVUo+kAKzH1OgDbfARCC5aoKFB2qBhQqbgjhbtIe6MgOw51tlAiSl7TNTn2I3wCEBCKhUKFVKXuNRZ7XbZCE1qlspW39pOOpFWhUWRGXI9CJPRdG0okeugKrA2O/LxE1Ip1GaOEFKDUKiQemYj0gH4OjyJzTIwxvQZYzYYY5z6BSWDskIrYxRis/w5caoabbyyTmu8T7gN9tCMhhQspOKhUCH1zHoAUTtaSRV4bDc0imwhL/pItkIrcAatEGmxKrIBTFXPEVLxUKiQusWqlTZA7SjrbXHEUy2SpQw2Qm4PNFlUhmDxETS5RjYxez53sqkIMgWdn72FkIqCU4pJXWFnUbWKyJAdEayH2jVSdqSwDunRSxKa+Asi4rz4J7G2kFZoYjAgrd7qRRZVl4gMeqcxWwGUtHWStk4rgFXOzDJbp6pzhJD6gaHvCSkRVlhE3bYbH18Z73FdAHqrPfkUqQ+o/iKkRNhEWqvcBvgAAiVi61GgkKqAIxVCSoyINAbNTlhIXUIqAQoVQgghRYPqL0IIIUWDQoUQQkjRoFAhhBBSNChUCCGEFA0KFUIIIUWDQoUQQkjR+P+VHEOwykuswAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "alphas = [1,0.85,0.7,0.55,0.4]\n",
    "De_cur = De\n",
    "for i in range(21):\n",
    "    if i % 5 == 0:\n",
    "        plot_cdf(100*mpcs, De_cur, label=f'$t={i}$', color = 'blue', alpha=alphas[i//5])\n",
    "    De_cur = law_of_motion.forward(De_cur)\n",
    "plt.xlabel('MPC (\\%)')\n",
    "plt.ylabel('Cumulative fraction')\n",
    "plt.legend(framealpha=0)\n",
    "plt.title(r'(b) $MPC_t$ distribution weighted by date-0 income, HA-one')\n",
    "plt.savefig('figures/figC1_b.pdf', format='pdf', transparent=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
